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A NOTE ON THE MOTION OF BUBBLES IN A 
HELE-SHAW CELL AND POROUS MEDIUM 


By SIR GEOFFREY TAYLOR and P. G. SAFFMAN 
(Cavendish Laboratory, Cambridge) 


[Received 12 August 1958. Revise received 13 November 1958] 


SUMMARY 

Exact solutions are presented for the steady motion of a symmetrical bubble 
through a parallel-sided channel in a Hele—Shaw cell containing a viscous liquid. 
The solutions also describe two-dimensional motion in a porous medium, since the 
two cases are mathematically analogous. It is shown that the solutions are not 
mathematically unique and a purely analytical criterion for specifying a particular 
solution is given. When the bubble is large, these particular solutions are such that 
the maximum width of the bubble is almost half that of the channel. Solutions are 
also presented for the motion of large asymmetric bubbles through the channel. 
The three-dimensional motion of bubbles in a porous medium of infinite extent is 
also considered briefly. 


1. Introduction 

THE motion of the interface between two immiscible viscous fluids in a 
Hele-Shaw cell, that is, in the region between two parallel, closely spaced 
flat plates, was recently described in a paper by Saffman and Taylor (1). 
In particular, they studied the propagation of a long finger of fluid through 
a parallel-sided channel in a Hele-Shaw apparatus containing a more 
viscous liquid. The boundary conditions at the interface depend on the 
interfacial tension as well as on the viscosities of the two fluids, but the 
exact nature of this dependence is unknown. To calculate the flow in 
the Hele-Shaw apparatus, however, only the change in pressure on passing 
through the interface and the thickness of the two layers of the more 
viscous fluid which are left behind adhering to the flat plates after the 
interface has passed are required. In (1) it was shown that the simplifying 
assumption that both of these are constant leads to an exact solution in 
closed form for the propagation of a semi-infinite finger. 

Actually, it was found that there exists a whole family of solutions 
corresponding to different values of the velocity of the finger and its 
asymptotic width. An experimental investigation showed that the width 
of fingers when the flow was not very slow was close to one-half that of 
the channel, and excellent agreement was found between the observed 
shape and the calculated one which was half the channel width. (At very 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959] 
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low flow rates, the observed shapes had a width greater than one-half that 
of the channel and did not conform with the calculated shapes having the 
same width. In this case it appears that the assumption that the pressure 
change across the interface is constant is not valid.) 

The primary purpose of the present note is to give the mathematical 
solution of the problem in which a bubble of finite size, i.e. a bubble 
bounded by a closed contour as opposed to the infinite ones considered 
in (1), moves steadily through a parallel-sided channel in a Hele-Shaw 
cell containing another viscous fluid, under the same simplifying assump- 
tions about the boundary conditions at the interface. Such bubbles some- 
times occurred in the experiments described in (1), although it was usually 
desirable to remove them, and the solution may not be without practical 
interest. However, it appears that the range of circumstances for which 
the solution is a reasonably close description of the actual motion is 
somewhat more restricted than for the long fingers. This is because, for 
reasons not yet fully understood, the physical conditions at the rear of 
the bubble, where the interface is retreating, may be different from those 
at the front where the interface is advancing, except in the case of very 
slow motion where very little fluid is left behind after the passage of the 
interface; but then the assumption of a constant pressure drop across the 
interface becomes less accurate. 

The analysis, therefore, is probably mainly of academic interest, but 
nevertheless may give a fair approximation to the shape of the front of 
the bubble where the assumptions are likely to be reasonable; this approxi- 
mation will be better the larger the bubble, and a long finger can indeed 
be regarded as the front of a very large bubble. Also, for the case of very 
slow motion and a bubble whose dimensions are small compared with the 
width of the channel, the pressure change across the interface due to 
surface tension and the interfacial curvature in the plane of the flat sides 
of the cell can be taken into account and the analysis should not be with- 
out physical significance. 

It is found that for bubbles of given area and for given conditions at 
infinity the solution is not unique, and that there exists a whole family 
of possible solutions. It is worth noting, however, that for the case in 
which the viscosity of the fluid inside the bubble and hydrostatic pressure 
gradients are negligible, a particular solution can be singled out by im- 
posing arbitrarily the condition that the product of the velocity and 
maximum width of the bubble is a minimum. For very large bubbles, 


whose front portion resembles a long finger, the width given in this way 
is half that of the channel, in agreement with the experiments, but a 
satisfactory physical explanation of this result has not been found. For 
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very small bubbles, the solution selected in this way is that bounded by 
a circular contour. 

The analysis of bubbles of finite size is restricted, for the sake of 
simplicity, to bubbles which are symmetrical about the centre line of the 
channel. Solutions will be presented, however, for the steady motion of 
very large asymmetric bubbles, that is, of long fingers similar to those 
considered in (1). The general asymmetric case seems to be capable of 
solution, but it did not appear to be of sufficient interest to justify the 
labour involved. 

The motion in a Hele-Shaw cell is mathematically analogous to two- 
dimensional! flow in a porous medium in which the motion is governed by 
Darey’s law. Simple closed solutions exist for the three-dimensional 
motion of bubbles in a porous medium of infinite extent, according to 
which the bubbles are ellipsoids of revolution. It is noted in section 7 
that the particular solution for which the bubble is a sphere is singled out 
by a minimum condition similar to that just mentioned above. 


2. Motion of a bubble in a Hele-Shaw cell 
In this section the analysis of the steady motion of a bubble in a parallel- 
sided channel in a Hele-Shaw cell is presented for the case in which the 


bubble is symmetrical about the centre line of the channel. We suppose 
first that the viscosity of the fluid inside the bubble is negligible and that 
pressure gradients due to gravity may be ignored (this will be so if the 
cell is horizontal). The mean velocity across the stratum of the viscous 
fluid through which the bubble is moving is given by 


b? op eh ap 
~ 1p ex ~ Be ~ ay’ 
(see, for example, Lamb (2) section 330), where z and y are coordinates 
in the plane of the plates bounding the cell, « and v are the components 
of mean velocity parallel to these axes, p is the pressure, y the viscosity, 
¢ = —(b® 12u)p is a velocity potential, % is the stream function (which 
exists by virtue of the equation of continuity), and the gap between the 
plates } is assumed sufficiently small for (1) to hold. (Equations (1) are 
satisfied by the two-dimensional flow in a porous medium of permeability 
b?/12.) 

Taking the walls of the channel as y = +1 and supposing the viscous 
fluid has unit velocity at infinity in front of and behind the bubble, ¢ and 
ys satisfy the conditions 


%’=+1 on y=+Il, o>x as ©£>+0 
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(edge effects which invalidate (1) within a distance of order 6 from the 
channel walls are neglected). 
To obtain boundary conditions satisfied by ¢ and y on the interface, 
we first assume that the pressure drop across the interface is constant. 
When the fluid inside is of negligible 
viscosity, the pressure inside is constant 
and hence 





¢ = const. = 0, say, (3) 


on the interface. The other assumption 
is that there is no film of viscous fluid 
left behind between the bubble and the 
plates bounding the cell, i.e. the ad- 
vancing interface expels all the viscous 
fluid in front of it, from which it follows 
that the velocity of the interface normal 
to itself is equal to the component of the 
mean velocity in the same direction. 
Expressed analytically, this is 








a 
on 8 és 

where n, s denote differentiation in the 
normal and tangential directions, re- 
spectively, and U is the velocity of the 
bubble. Integrating (4) and noting that 
the arbitrary constant of incvegration is 
zero because the bubble is supposed to be 
symmetrical about y = 0 gives, finally, 


Fic. 1. The motion of a bubble ,= Uy 

through a channel : (a) physical plane, ‘ 

(b) potential plane, (c) transformed 
potential plane. 


on ¢?= 0. (5) 
A discussion of these assumptions is 
given in (1). It is also shown there 
how, by a simple transformation, the solution for the case in which there 
is a film of constant thickness between the bubble and the plates may be 
obtained from that in which the films are of zero thickness. 

Now it follows from (1) that w = 4¢+i% is an analytic function of 
z = x+iy, and the solution of the problem is therefore effected by con- 
structing an analytic function of z which satisfies the boundary conditions 
expressed in (2) and (5). In Fig. 1(a) is shown the upper half of the 
physical plane; the maximum half-width of the bubble is denoted by A. 
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The potential plane, the boundaries of which are given by (2) and (5), 
is shown in Fig. 1 (6), corresponding points being marked with the same 
letter. Note that we must have UA < 1 in order that the potential and 
physical planes should be simple images of one another. The trans- 


formation tcos{ = tanh(47w)cot(47UA), 


where ¢ = £+in, transforms the potential plane into the semi-infinite 
strip » > 0, —4n < € < 4, shown in Fig. 1(c), corresponding points 
again being marked with the same letter. At the bubble interface, where 
¢ = 0 and » = 0, this transformation may be written as 


b= : tan-!(cos tan 47UA). 
7 


Writing z = w+2,, z,(f) = 2,(€,)+iy,(€, 7), we find that the boundary 
conditions (2) and (5) become in the ¢-plane, 


¥=9 on §=+h, (6) 


and ¥, = We 7 st ~ se im tanh-\(ie tan }rUA). (7) 


Here we have used the general formula 


; 2a 2b 
anh-"(a-+-ib) = 4tanh-! a eg 
tanh-"(a+ib) = }tanh (; io) +4itan- (i- aa 


It may now be verified directly that 


2= ioe 7 I tanh- ‘(ie tan UA) (8) 
7 


is regular in the strip and satisfies (6) and (7). The velocity potential and 
stream function are therefore given implicitly in terms of x and y by 
(expressing { in terms of w) 


. tanh-! canter”) (tanh? 4nw-+ tan? $7UA)!—tanh 4rv} |. 
n(4xUA) 
(9) 
Since the fluid inside the bubble is of negligible viscosity, the bubble will 
move faster than the fluid at infinity, ie. U > 1. Mathematically, this 
result arises as the condition that the transformation (9) be one-one (as 
may be seen without difficulty by considering the relation between x and 
¢ on the axis ¢ = 0). 
The interface is the image in the physical plane of ¢ = 0, —-UA< 4 < UA, 
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and has the equation (remembering that 4 = Uy and making use of the 


expression for the real part of tanh~1(a+-7b)) 


=~ tanh 1fsin?(4aU/A)—cos*(4aUA)tan*(4rU'y)}*. (10) 


The length BOD of the bubble is 


2U 
== 


4U 


— 
77 


l — . 
7 -tanh-{sin(47UA)} = L, say, 


and the area bounded by the interface, ie. by the closed contour given 
by (10), can be shown to be (see the Appendix) 
= c tanh MNtan?(drU A)! S, say. (11) 

The analytical solution (9) contains two parameters, U’ (the velocity of 
the bubble) and A (its maximum half-width), and it is clear that specifying 
the area S of the bubble provides only one relation between them. That 
is, if the area of the bubble and the velocity at infinity are specified, then 
there exists a whole family of solutions and possible bubble shapes corre- 
sponding to the various combinations of U and A which satisfy (11). The 
extreme shapes are A 0, for which UV 0, L x, and the bubble 
has zero thickness: and A 1, for which U 1 (since UA > 1) and the 
velocity everywhere is equal to the velocity at infinity. In the latter case, 
the bubbleis bounded by the channel walls and the lines x t+$L = +348. 
The shapes intermediate between these are ovals of length LZ and width 
2A, moving with velocities between U land U 00. 

If UA-> 1, keeping U fixed, the area of the bubble becomes very large 
and it is easily seen, changing the origin to the vertex of the bubble, that 
we recover the solution obtained in (1) for the motion of a long finger 
whose asymptotic width is A times that of the channel, namely, 

(1—A)log 3{1+-exp(—7vw)}. (12) 


3. The hypothesis of minimum UA 

Although the motion of a bubble in a channel does not possess, when 
formulated as above, a mathematically unique solution, we should expect 
it to be physically unique and that there is some mechanism, not accounted 
for in the analysis, to single out one particular solution. The search for 
some such effect has so far been unsuccessful, but several features of the 
analysis which distinguish a particular solution have been noticed. One 
of these is as follows. 
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Writing (11) in the form 


g — 16 MUA—-A), anh-Mftan(daUA)}, 
aw U*? 

it follows that S is a maximum for given UA and varying A when A = $UA, 
i.e. when U = 2. Alternatively, A is a minimum for given S and varying 
U’ when U = 2. Thus, if we make the hypothesis that the motion is such 
that UA is a minimum, then the bubble whose velocity is twice that at 
infinity is selected, and the corresponding value of A is given by (11) with 
U 2. 

The product UA of the velocity of the bubble and half its maximum 
width does not have a clear physical significance, although it may be 
identified intuitively with the rate at which fluid is pushed aside by the 
bubble. We have in fact been unable to place the hypothesis on a sound 
physical basis (attempts to relate it to a minimum rate of energy dissipa- 
tion fail), but it has the interesting feature that when the bubble is large 
and UA is close to 1, it gives a value for A close to 4, and this is the value 
observed experimentally with long fingers. 

It is perhaps worth mentioning, as a curiosity, that the shape for U = 2 
arises as the unique solution of a problem in the steady two-dimensional 
flow of an electric current in a uniform conducting medium. Suppose [see 
Fig. 1 (a)| that HF and AOE are two parallel electrodes at unit potential 
difference, separated by a medium of uniform conductivity, and that a 
perfect conductor of given area is laid on AE, the boundary of this con- 
ductor being BCDOB in the figure. If we now determine the shape of 
this conductor so that the amount of current entering it, i.e. flowing across 
the curvilinear boundary BCD, is a minimum, then it turns out (the 
details of the analysis are omitted for brevity) that the conductor has the 
same shape as the bubble of the same area for which U’ = 2. 

Another unique feature of the solution for a semi-infinite finger which 
is half the channel width may perhaps also be mentioned here. As has 
been pointed out by Zhuralev (3), the velocity potential due to the super- 
position of a uniform stream of inviscid fluid in a channel bounded by 
parallel walls y = +1 and a two-dimensional source of strength m placed 
midway between the walls at the origin of coordinates is 


w = qz+— log(2 sinh }7z), (13) 


a7 


where q is the velocity of the stream (a different coordinate system is 
used in (3) and the expression given there is accordingly different). 
It will be noticed that (12) and (13) would be closely similar if z and w 
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were interchanged in one of them. Now (12) may be written 
e-7W__ De-Mw/Al—Aomz/A1—-A) 4] re 


which can be solved explicitly for w when A = 4, giving 


w= 42’ re. log(2 sinh $72’), 
7 


where z’ = z—(log2)/z. This last expression is identical with (13) when 
m = 4q, except for a change of origin and a multiplicative constant. That 
is, the potential for a finger which is half the channel width is identical 
with that arising from the superposition of a uniform stream and source 
of appropriate strengths; and only the finger with A = } can be obtained 
in this simple way.+ (It can be shown that the potential for a finger for 
which A + } arises from a rather complicated system of sources, sinks, and 
dipoles superposed on a uniform stream.) 

Further, the dividing streamline of the flow due to a source in a uniform 
stream can be shown to have the equation 


ae lio ( sin y/A 
~— sin(1—A)zy/A)’ 


where 2A = 4m/(4q-+-m) is its asymptotic width; equation (13) can be 
interpreted as giving the flow past a fixed obstacle of this shape in the 
channel. Similarly, by taking axes fixed relative to the finger, we can 
deduce from (12) the flow past a fixed obstacle having the same shape as 
the finger, this shape having the equation (see 1) 


- 9(] — 
r= Lf i( + cos =) = =O tog(008 55) 


7 


When A = 4, the shapes of these obstacles are identical. 


4. The motion of small bubbles 
When A is small compared with 1 and U is finite, the dimensions of 


the bubble are small compared with the width of the channel, and (10) 


roduces to x2 
ee res | 14 
(U—1)? , i 


That is, small bubbles are ellipses with axis ratio U—1. The solution for 
a small bubble is obtained, in effect, by letting the width of the channel 
tend to infinity, so that (14) gives the shape of a bubble in an unbounded 


+ This result for the potential of the motion outside the finger remains true if the fluid 
inside the finger is not of negligible viscosity, the effect of gravity is taken into account 
and it is supposed that layers of fluid of constant thickness are left behind after the passage 
of the interface. However, the values of g and m will depend upon the particular 
circumstances. 
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Hele-Shaw cell. The corresponding complex potential is 


_ U—1, 2, nam, @ 
= (wt + UN +5 


The bubble of given size with a minimum value of UA is a circle of 
radius A moving with twice the velocity at infinity. Now the analysis 
has been based on the assumption that the pressure drop across the inter- 
face is constant. This pressure drop consists of a term 7'/R, where T' is 
the interfacial tension and R is the radius of curvature of the projection 
of the interface onto the plates bounding the cell, plus a term due to the 
curvature of the meniscus in a plane perpendicular to the bounding plates. 
A necessary condition for this assumption to be valid is that variations 
in T'/R should be small compared with changes in the pressure of the fluid 
surrounding the bubble, that is, small compared with 12uU R/b? (this 
condition may not be sufficient). When the bubble is small, this condition 
may not be satisfied, but it is interesting that the circular bubble remains 
an exact solution when the 7'/R term is taken into account, since R is 
then constant. 

This effect of surface tension tends to make the perimeter of the bubble 
as short as possible, and therefore to make the bubble circular. It is note- 
worthy that in this case this effect of surface tension and the condition 
of minimum UA give the same result. 

Visual observations of small bubbles which occurred during the experi- 
ments reported in (1) indicated that they were circular when sufficiently 
small. When larger, they were often pear-shaped with a rounded front 
and pointed back; some on the other hand seemed to be ovoid with the 
sharper end pointed in the direction of motion. These phenomena are 
probably due to the physical conditions at the retreating interface over 
the back of the bubble being different from those at the front; there is, 
however, no clear evidence on this matter. 


5. Bubbles of fluid with non-zero viscosity 
The previous analysis has been concerned entirely with the case in which 
the viscosity of the fluid inside the bubble is negligible and gravity forces 
are ignored. We shall now suppose that the fluid inside the bubble 
has a non-zero viscosity, that the z-axis is vertically upwards parallel to 
the channel walls, and that gravity forces are acting. The mean velocity 
across the stratum between the plates is now derived from a velocity 
potential be 
devs 12u 


where p is the density of the fluid and g is the acceleration of gravity. 


(p+ pg), 
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We shall denote quantities inside the bubble by the suffix 2 and quanti- 
ties outside by the suffix 1. The boundary conditions (2), (3), and (5) now 
become (still making the same simplifying assumptions about the physical 
conditions at the interface) 

wy +V on y= +1, d,>Vx as r++, (2)’ 


where the velocity at infinity is now denoted by V; also 


Pr a bp, qx = i b*p. gx (3)’ 
— 12u sien l2u 
wy bo Uy, (5)’ 


on the interface. 
It follows from (5)’ that the motion inside the bubble is given by 
¢, = Uz, yp, = Uy. (15) 
After some algebra we find that ® and ‘ defined by 


d, + ib, — [ [*z 


W—o+it ——*____, (16) 
V—U* 
~U Bb - 
where ? a C8> =. (Py—Po); (17) 
My 1 2p : 
satisfy (2), (3), and (5) with U replaced by 
U—U* = 
roe = U’, say. 


The expressions (9) and (10) with U’ replaced by U’ therefore give the 
relation between z and W and also the equation of the interface. In 
particular, it is to be noted that the family of possible shapes is independent 
of the physical properties of the fluid. 

The particular solution with U’ = 2 (and A close to } for large bubbles) 
is now obtained by minimizing U’’A for given bubble size. U’* is the velocity 
with which the fluid outside the bubble would move under the action of 
the actual pressure gradient inside the bubble, and the hypothesis of 
minimum lA has to be replaced by an hypothesis in which velocities are 
measured relative to |’*. The physical significance of this result is not 
clear. 


6. Penetration of an asymmetrical finger into a channel 

In this section we shall give solutions of the equations of motion which 
describe the steady motion of bubbles not symmetrical about the centre 
line of the channel for the case in which the bubble is very large and may 
be regarded as a semi-infinite finger. The motion of symmetrical fingers 
was considered in (1), and it was stated that there exists a singly infinite 
family of solutions for the problem of the penetration of a finger into a 
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channel filled with a more viscous fluid. This statement is in fact erroneous, 
and was made because it was not realized that further asymmetrical 
solutions existed. Actually, there exists a doubly intinite family of steady 
solutions, of which the symmetrical ones form a singly infinite sub-set, 
and it was thought worthwhile describing these here. 

It has already been noted that the question as to why long fingers 
formed experimentally are half the channel width, when their velocity is 
not too slow, remains unanswered, although several criteria of a purely 
mathematical nature have been discovered; the further question now 
arises as to why these fingers are apparently symmetrical, and a satis- 
factory answer has yet to be found. 

It will be supposed in the analysis that the viscosity of the fluid inside 
the finger is negligible and that gravity may be neglected. The solutions 
for the more general case may be obtained immediately by methods 
identical with those employed in section 5 above. (It is also assumed that 
no fluid is left behind adhering to the flat plates bounding the Hele-Shaw 
apparatus after the interface has passed by. Solutions for the case in 
which films of constant thickness are left behind may be derived in the 
manner described in (1).) We employ the same notation here as in section 
2, except that the origin of coordinates will no longer be the centre of the 
bubble. 


The velocity potential satisfies 


d=—0 (3) 


on the interface, by virtue of the assumption that the pressure drop 


across it is constant. For the stream function, the boundary condition (4) 


ob U(y—Yo) (18) 
where U is the velocity of propagation of the finger and y, is a constant 
whose value is as yet unspecified. (y, may be taken positive without loss 
of generality, solutions with y, negative being mirror images in y = 0 of 
those with y, positive.) Taking unit velocity at infinity ahead of the finger, 
we have the further conditions ¢ ~ 7, J ~ y as x—> +. 


now integrates to 


The symmetrical solutions are obtained by putting y, = 0. It was shown 
in (1) that the complex potential for the steady motion of a symmetric 
semi-infinite finger is given by equation (12), where A is the ratio of the 
asymptotic width of the finger to that of the channel and the velocity of 
propagation U = 1/A. The equation of the interface (i.e. the image of 
¢ = 0) may be written in the form 


ore. S7o3 
r= 21 Oo hog(oos jnl'y), (19) 
7 
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the origin of coordinates being taken at the nose of the finger. The poten- 
tial plane (see Fig. 6 of (1)) is the semi-infinite strip ¢ > 0, —1 < %< 1; 
the line ¢ = 0, —1 < %& < 1 is the image of the interface, and the lines 
d > 0, 4 = +1 of the channel walls. 

The method by which the symmetrical solutions are obtained may be 
extended with little difficulty to deal with the asymmetrical case in which 
the boundary condition (18) applies. We have that x and y are conjugate 
harmonic functions of ¢ and % such that y= +1 on Y= +1, d> 0; 
y~d as d6>+oo; and y= U-Y+y, on d= 0, —1 << 1. Take 
therefore 

y = y+ s A,,e-"7¢ sin nmb+ s B,, e-"*™4+ sin 4na(p+1), 
n=1 n=1 
thereby satisfying the boundary conditions on ¢ = +1 and at ¢ = o. 
The boundary condition on ¢ = 0 is satisfied if 


U-b+y, = b+ 3 A, sinnw+ > B,sin}nr(p+1) (—1 << 1). 
1 1 


Calculating A, and B,, by the usual method of Fourier series, we find that 


2 . 2y 
A, = -(] — T -1), B,, = fl —(—1)"}. 
ated nt 
The series for yd r- Ly may be summed, giving finally 
2 . 2y 1 +ie-t7 
2 P . | g—mw) 4 “70 ie Be ‘ 


The interface is the image of 6 = 0 and has the equation (remembering 
that & U(y—Yp)) 


) » 
x = —(1—U-)log cos 4nU (y—yy) + “Vlog tan{}a+ knU(y—yp)}. (21) 
1 7 


Now it may be verified that the transformation (20) is one-one and free 
from singularities if ’ > 1 and y,+U-! < 1. The former of these condi- 
tions is clear physically because the fluid in the finger is less viscous than 
that in the channel. The latter condition arises from the fact that all the 
streamlines must originate on the interface (simply by continuity, because 
the velocity of the fluid at x - is zero since conditions must be 
uniform along the straight sides of the finger; hence ¢ = 0 for all y at 
x —oo). Hence, the value of y at the point where the streamline y = | 
meets the interface must be less than one, and by (18) this value is y,+ U-?. 
Y, is in fact the value of y at the point where the streamline ¢ = 0 meets 
the interface. 

As regards the shape of the interface given by (21), we note that as 
y>yot+U4,2~ (1—U—y,)log(y,+ U-!—y) > —oo; and as y> y,—U-, 
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x~ (l—U-'+y,)log(y—U-!+-y,) > —o. The interface therefore has the 
two asymptotes y = y,+U-!, y = y,—U-—, and lies between them, being 
nearer the wall y = 1 than the wall y= —1. The asymptotic width of 
the finger is 2U-! = 2A, say. The interface for the values U = 2, ys = } 
is sketched in Fig. 2. 


y=, ws =1 











y=, w= -1 
Fic. 2. An asymmetric finger penetrating a channel. The interface is 
calculated from (21) with U = 2, yp = }. 


If y,+U-! = 1, the interface intersects the wall y = 1 at right angles. 
In this case, the interface is half of a symmetrical finger between the walls 
y = —l andy = 3, as is clear from the fact that (21) may then be written 
in the form 


ta =x 2 — U-)log cos }nU(y—1). 


7 


There is no feature of the analysis which specifies the values of A (or U) 
and y,, and therefore there is a doubly infinite family of mathematical 
solutions for the penetration of a finger into a channel containing a more 
viscous fluid. It remains to explain why the fingers obtained in the 
experiments described in (1) appear to conform with the shape calculated 
from A = 3, ¥% = 9. 


7. Bubbles in three dimensions 

Provided that the same simplifying assumptions are made about the 
physical conditions at the interface, the motion in three dimensions of a 
bubble in a uniform unbounded porous medium, filled with viscous fluid 
moving uniformly at infinity, and in which the motion is governed by 
Darcy’s law, possesses simple solutions in which the interface is an ellipsoid 
of revolution. These solutions are briefly referred to in an article by 
Polubarinova-Kochina and Falkovich (4), but since they are apparently 
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not well-known, it was thought worthwhile to present the results in a 
more general form than hitherto, and to point out the lack of uniqueness. 

If V denotes the velocity at infinity, U the velocity of the bubble, g the 
acceleration due to gravity, and 
. k M211) 
U® = Si... Pi) he e 

My J 

where k is the permeability of the medium and the suffixes 1 and 2 refer 
to the fluid outside and inside the bubble, respectively, then the boundary 
conditions that the pressure and normal velocity are continuous across 
the interface can be satisfied by taking the bubble as any ellipsoid of 
revolution with axis parallel to U—U*. The relation between U and V is 
3 


jog. Pr —e\(U—U*) -(V—U*) (prolate ellipsoid), 
f e~ 


(1—e?)' sin-'e)(U— U*) e3(V—U*) (oblate ellipsoid), 


(¢ 
where ¢ is the eccentricity of a meridian section. 

The solution for a bubble of given volume is again mathematically not 
unique, but a particular solution is singled out by the hypothesis that 
the product of the equatorial area A and the velocity measured relative 
to U* should be a minimum for a bubble of given size and for a given 
velocity at infinity relative to U*, i.e. that 

A|U—U*|/ V—U* 
isa minimum. The bubble shape is then spherical and U— U* = 3(V—U*); 
this solution remains exact if interfacial stress effects are present whose 
effect is to produce a pressure drop proportional to the interfacial curva- 
ture. There does not appear to be at present any experimental evidence 
on the motion of bubbles in porous media. 


APPENDIX 
To calculate the area bounded by the contour given by equation (10), we note 
that on the interface (where @ 0 and us Uy), 
xr r—o r and Wy y us . y(U 1). 


Hence, we have from (8) 


4U-1 , , 
xr—(U—l)y j7 tanh ie" tan(}7UA)} 


4 


4 U—1 XS (—1)tel@*€ tan?*+1(47UD) 
we U hw 2n+ 1 


77 
n=O 


on the interface (remembering that 7 0 there). The equation of the interface may 
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therefore be expressed parametically by 





, 


2n+1 


4U-1 S (— yn Sin(2n-+ DE tan®™+(JarVA) 
t= -- - - 
a U 
0 





+ 
2n+1 


x 
4 cos(2n + 1)€ tan?"+1(42rUA) 
y= (—1)" 
0 
and the complete contour corresponds with —a < € < 7. These series are uniformly 
convergent if tan(}7UA) < 1, which is the case if the bubble dimensions are finite. 
We have now 


I-1S tan*"*2(4arUA) 
0 


16 U-—1 ; — 
te tanh~!(tan }7UVA)’, 


by virtue of the relations 


7 

| cosm€cosn€ dé = 0 if m 
=e . 

7 ifm 
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MOTIONS WITH INITIAL GRADIENT 


By Sr. I. GHEORGHITZA (University of Bucharest) 
[Received 22 May 1958] 


SUMMARY 


This paper is concerned with the motion of incompressible fluids in a porous 
medium taking place only when the modulus of the gradient of the dynamic pressure 
surpasses a given value. The equation (12) which is satisfied by the function ¢ defined 
by (1) and (3) is deduced. 


A general approximate method of solution is described, and some particular exact 
solutions are obtained. 


1. Introduction 


A DETAILED study of the motions in porous media shows that Darcy’s 
law, or the linear law, expressed by the relation 


v k grad #3 ) —k grad H, say, (1) 
in the case of isotropic media, is valid only over a limited range of veloci- 
ties. In particular for very small velocities the fluid is moving more slowly 
than is required by (1). Consequently, a limiting velocity may exist 
under which (1) ceases to be valid. We must mention that, at the same 
time as Darcy, Schmidt (see 1), by studying the motion through animal 
membranes, found that fluids can pass through the membranes only after 
the pressure gradient exceeds a certain value; until this value is reached 
the membranes are practically impervious; but an animal membrane is a 
particular porous medium. Recently, Swedish and other scientists found 
the existence of a limit gradient at which fluid motion in the ground stops 
(2). This limit gradient will be named initial gradient; motion occurs 
only when this initial gradient is exceeded. 

These remarks show that it is not without interest—theoretical and 
practical—to study the motions in porous media with initial gradient, as 
defined above. The existence of the initial gradient shows that in certain 
conditions account must also be taken of the rheological aspect of the 
motion. 

In this work we shall derive the equations which describe the motion 
in porous media, illustrated by some simple examples, and shall apply a 
classical approximate method. A result worthy of mention is that, in the 
plane case at least, the discharge below a dam on a homogeneous unlimited 
porous medium is no longer infinite. 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959] 
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2. The equations of motion 


Mathematically, the motions with initial gradient can be expressed as 
follows: 


v= 0 for \gradH| < K* 


K* grad H rs 
grad H for grad H| > K 


where the critical value of \grad H| is denoted by K*, that is, the magni- 
tude of the initial gradient. These equations show that in the domain 
filled by the fluid there may exist regions in which there is motion, and 
regions in which there is no motion. We shall confine ourselves to the 
case of homogeneous and isotropic media and to incompressible fluids. 
Introducing the velocity ‘potential’ ¢ and the constant K, defined by the 


relations ¢ = —kH, K = kK*, (3) 


2 
Vv —K( grad H— , (?) 


(2) may be written 


v= 0 for \grad¢d| < K, (4) 
v= (1- = — }erad ¢ for \grad¢| > K; (5) 


when grad 4) is sufficiently large, (5) tends to (1). The relations (4) and 
(5) will be the basis of our study. The relation (5) is valid for steady 
motion, but it can be generalized for unsteady motion to the form 

. CV 4 

é ne (1 niger e—v: (6) 
for K = 0 we obtain again the general vectorial equation of the motion in 
porous media (see for example 3), while for K + 0 and steady motion we 
obtain the equation (5). 


mg et 


Taking, as in the classical case, 
k év 
V=v+4— —, (7) 
mg ct 
the equation (6) becomes 


a : 
V= (1 ~ aap) red. (8) 


There is no reason for introducing the vector V defined by (7) in the 
regions where \grad¢| < K; in these regions V and v must be zero. But 
in the regions where \grad¢| > K, we shall have from (7) and (8), and 
assuming that ¢ has been determined, 


a 


Vv Vee matt 4 Fema { onan — 


gradgdt. (9) 
: U 
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3. The equation for ¢ 

It is known that in the case of Darcy motions ¢ is a harmonic function 
in the domain occupied by the fluid. By using (5) or (8) and the con- 
tinuity equation, we can derive an equation for ¢ in the regions where 
the motion takes place; further the equation 


div Vv 0 (10) 
implies that div V = 0. From (5) and (10) we obtain for ¢ the equation 
K 


grad d 


‘ K . 
|V*¢ T grad ¢. grad i" — 77) = J (1 1) 


(12) 


‘i v4 ; K (v96— grad d.grad paced -\, 
grad ¢@ \ 2 grad |? 


Putting K = 0 we obtain V7¢ = 0. 


We need now an equation which gives us ¢ for |grad¢| < K. The most 
natural hypothesis would be one of the following: in the regions in which 
grad¢| < K, either ¢ satisfies the equation (11), or it is a harmonic 
function, as in the case of Darcy motions. Here we shall use the first 
hypothesis, because the use of the second one presents very great diffi- 
culties. In particular, we would have to impose some boundary conditions 
connecting the functions é and ¢ of the two subdomains: D,, the domain 


in which there is motion, and the domain D, in which there is no motion. 


First, we would impose the continuity of the pressure, which is expressed 


by db? = g® onS, (13) 
S being the common frontier of the two domains. Secondly on the 
boundary separating the two domains we must have the limiting value of 
grad d™ |, i.e. 

grad 6" K onS. (14) 
The boundary of the domain in which the velocity vanishes would then 
behave for the fluid in motion like the surface of an impervious body for 
a viscous fluid. 

However, with the first hypothesis mentioned above, (13) is satisfied 
automatically, and (14) is the only condition requiring to be satisfied on 
the boundary of the domain in which there is motion. 

There are also conditions to be satisfied on the boundary of the domain 
filled by the fluid in the porous medium. But these are the ordinary 
boundary conditions relative to the function 4, which are found in any 
special treatise on underground hydrodynamics. 

In this way the problem is reduced to the integration of equation (12) 
with certain boundary conditions, the boundary of the domain where we 
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have motion being determined by the condition (14), and the velocity 
being given by (5) once ¢ is known. 
4. An approximate method 
To solve particular problems a classical approximate method is some- 
times useful. Put wear . 
¢ = $)+K¢d,+K%db.+..., (15) 
where K is a small parameter having the dimensions of a velocity. Now 
introduce a dimensionless parameter A, defined by the relation 
K, = KV;*, (16) 
VY, being the modulus of a characteristic velocity which is actually realized 
in the domain D,. We shall always have AK, < 1, and we may write 
, 2 - 
$ = b,0+ Kodo: + Ko do2+-5 (17) 
do; (j = 0, 1,2,...) denote the functions corresponding to the expansion 
in terms of Ky. It is evident that 
$; = Veiga; (j = 0,1,32,...). (18) 


The equations which determine ¢5, ¢,, ¢5,... are immediately deduced by 

the substitution of the expansion (15) in the equations (12). We find 
V24, = 0, (19) 
v4, grad d,.grad grad ¢, ” 


2 grad ¢,|* 





; (20) 
—_____— | grad d, grad dg. grad grad ¢, |? + 
2 grad ¢,|* : s , 

- grad ¢,.grad grad ¢, ?+-2 grad ¢, .grad(grad ¢,.grad ¢,)— 

—2 grad ¢,|~'(grad dy. grad ¢,)(grad ¢,.grad grad ¢,\?)| (21) 
and so on. After solving these equations, subject to their boundary condi- 
tions, the boundary of D, is found from the relation 

gradd| = K, (14’) 
that is 
grad ¢, +A grad dy ~'grad dy. grad ¢,+4K*(\grad ¢, ?|\grad ¢y|-!+ 


» 


+ 2'grad dy —'grad $9. grad ¢,—(grad $y. grad ¢,)*|grad ¢y|~?+ 
+... = K. 


Then, from (5), the velocity is 


. rad dy \ , 
v = grad¢e4 K(grad 4, — Bras fe) 


” d dy. grad 
+K *(erad $2+ aire grad $4) Foes 
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The first boundary-value problem determines a harmonic function and 
the others determine functions which satisfy the Poisson type equation. 


5. Examples 


We shall now consider some important examples of motions for which 
exact solutions are obtainable. 


(a) One-dimensional motion in the x-direction, of a layer stretching 
from x = 0 to x = /, with p, and p,, denoting the pressures at the ex- 
tremities. We obtain for the velocity 


u = K(py—py,)(pgl)-'—K (23) 


(6) A plane motion with spherical symmetry, denoting by p, and py, 
the pressures for r = r; and r = r,, (taking ry; > r;, py > py, to fix our 
ideas) and by q the fluid discharge for unit thickness of layer; from the 
equation 


q _ dd _K. 
2rr dr 
we obtain at once 
1 
q 2n[ k( p,- Pr) pg) i. K(rn—n)]( in) . (24) 
I/ 


(c) For three-dimensional motion with spherical symmetry, and using 
similar definitions to those in example (6) we have 


Q = 4a[k(py—pPr (pg) —K(ry—rp (rp ri) 3; (25) 


here Q represents the total fluid discharge. 

From the formulae (23)-(25) it follows that the discharge diminishes for 
the motions with initial gradient, as was to be expected. 

(d) Let us take now a porous medium with the following boundary: 
a<—R,y=0;r>R,y=0;-R<ax< R,y = ,(R*—2%), the y-axis 
being directed vertically downwards and in the medium y > 0. We con- 
sider the rectilinear parts of the boundary as supply surfaces and the 
curvilinear one as an impervious surface, so that 


“a (0 fory=0,2r>R 
\|-—kH fory=0,7%<—R 
ch — , 
re 0 ony=,(R*—2x*?) (—-R<2a< R). (26) 


It is easily seen that the exact solution of the problem is given by the 
harmonic function 

kH®@ 
— ’ 


_ 
77 


d (27) 
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where @ is the angle made by the z-axis with the radius vector to a 
given point. The velocity will then have the expression 

ae -(4 —K}e, (28) 

mr 

where € is a unit vector normal to the radius vector. We find that the 
boundary of the domain in which there is no motion is a semicircle with 
centre the origin, and of radius 
Sk 


* 
. ak’ 


(29) 
that is, v = 0 forr > R*. 


On calculating the discharge of fluid from upstream to downstream, we 


find 1 
= <eIn(e)—K(R*—B), 
7 R 
kH kH » 
or Q — a {in ( ae) 1} + kK R, (30) 


so that we obtain a finite discharge; if K + 0, then Q > 00, as for Darcy 
motions. 


(e) Finally, let us consider the following unsteady one-dimensional 
motion: 
* 


t< 0, pO)= pt, P(L)=py PT > Pu 
t>0, p%)=p, PL)=py Pr>Pu Pr #PT 
As p = ax+b (a, b are constants), it follows that we have for p 


(31) 


(<0, p= pit+(Pu—Ppel; 
t> 0, p= Pr+(Pu—pprL-. 
The equation of motion will then be 


dv mg, ,m re _pgKL 


a ‘k “pp L k(py— Pur) 


of which the solution is 


an 0, (32) 





— k KL 
v=v mot 4 Pr Pur te _ PY 1—e-motk) 33 
. L pg ‘(Pi—Pu) : oe 


. at pal pg KL 
where v, is Cg SE tenes fee ge . 34 
’ " pgL k(p?—Pu) (34) 





For t > «, we have 
ime = # Pree KL 
pg =OCUL kK(py—Py) 





as should be the case. 
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(Armament Research and Development Establishment, Fort Halstead, 
Sevenoaks, Kent) 


[Received 15 May 1958] 


SUMMARY 
The behaviour near the sonic plane of an axially symmetric sonic jet of gas, 
emerging into a region of lower pressure, is studied. The form of the singularity at 
the sonic plane is found, and the lowest terms of the series expansions of the variables 
of state near to the plane are obtained. In particular the corresponding solutions for 
this case and for the two-dimensional problem are compared, and it is shown that 
there are important differences between them. 


1. Introduction 

Ir a gas emerges as a jet from a plane circular orifice into a region of 
lower pressure, in such a way that the velocity in the plane of the orifice 
is uniform and supersonic, and so that the streamlines are everywhere 
normal to this plane, then the jet will expand, and the velocity will increase. 
The behaviour of the jet near the orifice has been studied by Johannesen 
and Meyer (1). An expansion wave arises from the edge of the orifice 
which interacts with itself, and is reflected from the jet boundary, giving 
rise to a series of interactions and reflections (see Fig. 1). 


Orifice 
—~ 
M>1—~> “SSE 
one ss P 
SVT OS 
>. eS 
gp > == » 
Saat Sr 
~ 
4 ~ ~ 
Axis ee a 








Fic. 1. Initial characteristics for flow from a slightly supersonic orifice. 


If the velocity at the orifice is sonic the situation is more complex. 
The characteristics which form the upstream limit of the expansion fan 
at each point of the edge, all lie in the plane of the orifice. Thus instead 
of a point intersection, they intersect at each point, and the region of 
interaction begins at the orifice instead of downstream, as in the former 
case (Fig. 2). It is clear, therefore, that in this case a singularity must 

+ Now at Sir John Cass College, Jewry St., London E.C, 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959) 
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occur in the plane of the orifice; the study of this singularity forms the 
content of this paper. 


Orifice 








Fria. 2. Initial characteristics for flow from a sonic orifice. 


The singularity is not of course peculiar to the axially symmetric 
problem. The equivalent two-dimensional problem was studied in (2). The 
singularity in this case was representable in terms of the Weierstrassian 
elliptic function. It is of interest to compare the two solutions, and this 
has been done in this paper. 

The equations for the subsequent terms in the expansion will be non- 
homogeneous ordinary linear differential equations of the second order. 
Thus in principle a numerical solution to any order of accuracy can be 
obtained from this solution. 


2. The equations of motion 
The equations of motion of a homentropic perfect gas in axially sym- 
metric steady flow can be expressed in the form 


Cue 
= 0 (2.1) 
Cr Oz 

, , é , ., ov .v 
u — {log(p/p,)}+v — {log(p/p,.)} +—+—+-=9, (2.2) 

2 or 2 Or ? 

2 2 y+! , 
usr + a* = —— a, (2.3) 
a y= 

(ajay)? = (p/py)’-, (2.4) 


where conventional notation has been used, and where the starred values 
are those reached on the sonic line. 

Consider the flow in a meridian plane. It is convenient to transform 
to coordinates with origin at the upper lip of the orifice. If the radius of 
the orifice is r, let r3¢2 = 284-(r.—r)P, (2.5) 
and {r,—ritand = z, (2.6) 


as shown in Fig. 3. Let Ua, be the component of velocity in the direction 
of increasing €, and Va, be the component in the direction of increasing ¢. 
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Let A = a/a, and » = log{p/p,}. In terms of the new variables equations 
(2.1) to (2.4) become 
gf v2 8 








7 
of | ch > (2.7) 
> OF , oU av & = ; 
Lp OG ne OW A Se —U cos¢} = 0, 
é cf Va6 oe bag +9 saad ane mers 
' (2.8) 
: = 2 y+1 
Ut4 y24 —— A? — 2, 2.9) 
y—l1 y—1 
log A = sa (2.10) 
| 
| = 
% 
7 
1 
! 
Axis : 








Fic. 3. The coordinate systems. 


The interesting region is that for which ¢ is small. The obvious ap- 
proach, therefore, is to expand the dependent variables in powers of ¢. 
The symmetries of the problem show that the appropriate expansions are: 

A = 1+A,(£)¢?+A,(E)b*+ O(¢), 
V = 1+V,(£)g?+Vy(E)d*+- O16), 
U = $+ UE)? + UlE)d? + O14"), 
7 = nal€)P? + nalE)b*+ O19"). 

The boundary conditions are that at the edge of the orifice the flow 
must be locally a Prandtl-Meyer expansion, and, since the coordinate 
system is biased by being rooted on one lip of the orifice, that the solution 
must be invariant to a change to a coordinate system based on the other 
lip. 

The first condition is most conveniently expressed by 

A=V at €=0, forall ¢. (2.11) 

The second condition consists of four equations of the form 


n(E1; $1) = (Eo; be), (2.12 


where the suffixes refer to opposite edges of the orifice. Only one of these 
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equations is actually required as a boundary condition, the others, with 
the equations of motion, being satisfied identically. These two coordinate 
systems are related by the equations 
& 4+ £?—4£, cos dy, 
£, sind 

and tang, = —"1 1 _.. 
2—€,cos¢, 
Expanding these in series for small ¢,, they become 

5 X{1 T Es {dy X}? +-O(¢4), 
and d, = €,{¢,/ X}[1—3{X2+3¢, X+ +-O(¢$), 


where X = 2—£,. Substitution into (2.12) gives 


X?yo(E,)675 = EF no(X)b3+ O(¢}). (2.13) 


3. The first approximation 

If the series expansions are substituted into equations (2.7) to (2.10), 
and the coefficients of the lowest powers of ¢ are equated to zero, only 
two independent equations are derived, namely 


2(.+53)+1 0, 


and A, 


It will easily be seen that, in general, from the set of four equations 
obtained from the coefficients of successive powers of ¢, the functions 
with highest suffix can be eliminated to give a non-linear relation between 
the functions with lower suffixes. 

Thus the next powers of ¢ give 


dn r 
ae 


and 


Eliminating ,, Vy, and A, from (3.4), (3.5), and (3.6) we find 


1 2noV, gts _ su, 


< > dé “ 


.dn,» 
g dé 








BEHAVIOUR OF AN AXIALLY SYMMETRIC SONIC JET 291 
If this is combined with (3.1), (3.2), and (3.3) an equation for », is derived, 


namely, 


24°» - dn. 2g + 6(y-+- = 
C de — (s+ 5 )eR (6+ = 2)n + (y+ 1)n3 0, (3.8) 


and this can be simplified by the substitution 


(y+1)n, = —&F(€), 


to the form (1—é)F"’—F’ = 6(1—&)F?. (3.9) 
Using (2.11) and (2.13), the boundary conditions become 

&—> +0, &F(é) > 1, (3.10) 
and F(é) = F(2—&). (3.11) 


As pointed out in the introduction, a problem equivalent to the two- 
dimensional problem was solved in (2). The analysis is very similar to 
the above, but instead of equation (3.9) the resulting equation is 

F" = 6F?, (3.12) 
though the boundary conditions are identical. Thus in the usual notation 
for the Weierstrassian elliptic function the solution is 

F(é) = 29 (a€), (3.13) 


where « is chosen to give the function a period of 2. 


4. Solution of the fundamental equation 

It is easily established that equation (3.9) has no solution in known 
functions. Near € = 1, almost all solutions of (3.9) behave like 

F() = k1—€)*+0(1), 
and so F(é) has a double pole at this point. In this problem, however, 
the only possible singularities which could arise physically are at € = 0 
and € = 2. Thus since F(€) has no singularity at £ = 1, k = 0, and since 
further F(é) is an even function of (1—£), the boundary condition (3.11) 
can be written F’'=0 at €=1. (4.1) 
If in (3.9) we write 
l—ar=€, and Af(x) = F(é) 


and choose « and £ to satisfy 
a8 = I, 
then the equation becomes 
d gs 
4 af’ a) = Gap? (4.2) 
d 


If further yb = log x, y = xf (x), and A = i’ 
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then (4.2) reduces to the first order equation 


A da _ 4A - 4y —6y? om ©. 
dy : 


The boundary conditions are that 
A=0 when y=0, 
and is a regular function of y near the origin, and secondly that 


(bo—vb)*y¥>1 as fy, = FlogP. 








Fic. 4. Integral trajectories for small y. 


\s this implies that 
(bo—P)PA>2 as pop, 
an equivalent condition is that 
A= 2y!i+o(y!) as yoo. 
The finite singular points of equation (4.3) are at the points 
y == 6, A 0, 
and at y = 4, A = @, 


It is easy to see (Fig. 4) that the first point is an unstable node, while the 
second is a saddle point. Only one of the integral curves through (0, 0) 
is regular in the region around the singular point. This is the one which 


behaves like r 


2y 
in the neighbourhood of the origin. Since, physically, our solution must 
be regular at this point, this is the required solution. 
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In the neighbourhood of y = 0 therefore a series expansion of the form 


— ll+ 5 
A= 2y{l+ Yay (4.4) 
is valid. 
Substituting into the equation, and equating to zero separately the 


coefticients of the various powers of y, the following recurrence relations 
are obtained: 


(4.5) 


ar 
and if r >: 2 : > (4.6) 
s=1 
so that the coefficients alternate in sign. Now if r > 2 


> 


2r >r+22>Fr, 


oot | 
and so > a, |< < > |a,|\a 
s=1 


x 
Consider the series L,(y) = 1+ > A,,y’, 
r=1 


and L*(y) = 1+ > Afy’, 
r=1 


where Ay, = Aj = 4, 


and for r > 2 


r—-1 r—1 
A i j 2 Aas Aer» A? = 2 At Are 


Then for all r Ase < \a,| < A®. 
If R, is the radius of convergence of L,, R* that of L*; and if R is the 
radius of convergence of the series contained in (4.4), then 
R* 4 R << R,. 
But 
ee r- 


J 1 
L2(y) 1+ > 24 *r > A,,A oro} 


1 
= 1+24,,y+4 Any = 4L,(y)—2a,y—3. 
r=2 
Thus L,(y) = 2—{1—2a, y}}, 


and R, = as 


° 2a, 
It can therefore be established that 
} => Ra, > }, 


and so the series contained in (4.4) converges uniformly within a circle 
whose radius is at least }. 
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To study the other end of the range it is convenient to put 
A=2y'T and yf?=1. 


,aT 


Then (4.3) gives Cl i 
ag 


iT 








Fic. 5. Integral trajectories for small Z - 


The singularities of equation (4.7) are at the points 
C 0, r +1, 


and C +3, r U. 


The first two singularities are unstable nodes, and the latter are saddle 
points. Fig. 5 shows the form of the integral curves. Near ¢ = 0, the 
solution will behave like the solution of 

av 


i 6([—1)+4Z = 0, (4.8) 


and it is at once obvious that the integer coefficient of the second term 


will be troublesome. Put 
P() l+ of T oD C* oF; CF - oD ¢*+ 050° +y(¢), 
and choose the coefficients 9b,, gbs,..., 9/; 80 that when this expression is 
substituted into (4.7), the resulting terms in the first five powers of ¢ 
vanish. Then 
6y+4,b, [®+-0(f*) = 0, (4.9) 
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and so y = —4 yb, C% log [+ 9b, *+-0(C°). (4.10) 


Thus the solution of (4.7), valid for small f, is 


P(t) = > > .b, L%**{log ¢}s. 
r=0 s=0 
If in this expansion we have 
o% = 1, 


and yoo +49); = 9, 


then all the other coefficients are determined by the recurrence relations 
obtained by substituting (4.11) into (4.7), with the exception of 9b,, and 
there is no recurrence relation which determines this. Thus 4b, appears 
as an arbitrary constant, and the later coefficients are functions of it. 
In fact it may easily be shown that the coefficient ,b, is a polynomial in 
o¥¢ of order [r/6], where [2] means the integral part of x. 


5. The integration and the results 

Since it has been shown that the series (4.4) is uniformly convergent 
within a finite region of the origin, within this region (4.4) may be used 
to give numerical values not only for A, but also for its derivatives of all 
orders. 

The trajectory chosen avoids the other singularities in the finite part of 
the plane, and so equation (4.3) may be integrated numerically to as large 
a value of y as is required. 

Unfortunately it does not seem possible, by elementary processes, to 
establish the convergence of the series (4.11). An intuitive argument for 
its validity can however be found. If the numerical solution of (4.3) 
for large enough y is equated with the solution (4.11) of equation (4.7) 
for the corresponding value of ¢, the unknown yb, can be found. If then, 
as does in fact happen, different values of y give the same value of 9b, the 
process is justified. 

Thus by the above process the function 7,(€) was calculated, and is 
plotted in Fig. 6. Also in this figure for comparison is plotted the corre- 
sponding two-dimensional solution. It will be seen at once that there is 
a fundamental difference between the two curves, as the two-dimensional 


solution not only passes through the origin but does so with zero slope, 
while the axially symmetric solution has a finite slope at the origin. Also 
the value of »,(1) is greater than 2, whereas the corresponding value of 
the two-dimensional solution is less than 1-5. 

After the non-linear equation for the function »,(£), the equation for 
n,(€) and the higher functions will be non-homogeneous linear equations. 
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The equation, however, and its boundary conditions are rather tedious 


to set out, and so are omitted. 


Axially symmetric solution 


Two-dimensional! solution 


1)7), 


" 
w 


(y+ 
/ 








Fic. 6. A comparison of the variation of —(y+ l)log(p/p,) with é, in 
the two-dimensional and axially symmetric flows. 


6. Conclusions 

An approximate solution has been computed for the flow in an axially 
symmetric sonic jet near its sonic plane. In particular the singularity in 
the sonic plane has been examined. The solution obtained here could be 
used as the starting conditions for a characteristic net, which would 
determine the flow through the sonic jet. 
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SUMMARY 


A solution of Tricomi’s equation is obtained for the flow past a thin, doubly 


symmetric body placed at zero incidence in a high subsonic stream in which sonic 


velocity is attained along a segment of the body. This fiow is the compressible 
analogue of the Riabouchinsky model for incompressible fluids. The singularity in 
the hodograph plane corresponding to the point at infinity in the physical plane is 
essentially different from that which occurs in other similar problems. The boundary 
value problem is of mixed type and this is shown to lead to a pair of dual integral 

juations for which the solution is obtained. Numerical results are given which 
specify the dimensions of the body corresponding to a range of incident Mach 
numbers. By symmetry the total drag on the body is zero. 


1. Introduction 


THE study of problems in transonic flow in two dimensions is frequently 
based on solutions of Tricomi’s equation 


> u 
cr" cry 


u— = 0, (1) 


C ua or 
where u and v are related to the velocity components and y is the coordi- 
nate measured at right angles to the main stream. This equation is an 
approximation to the full hodograph equation of a perfect gas and is valid 
for small perturbations about a uniform flow at a Mach number of unity. 
In recent years a number of solutions have been found to represent the 
flow past a symmetric body whose axis of symmetry is parallel to the 
direction of the main stream when the flow at infinity is sonic or slightly 
subsonic. These bodies usually have straight sides but they may also 
consist of portions on which the magnitude of the velocity instead of its 
direction is constant. Such problems may be formulated precisely in the 
hodograph plane and the corresponding solutions obtained. For wedge- 
shaped bodies the sonic point is taken at the shoulder of the wedge. 
However, the continuation of the flow beyond this point presents diffi- 
culties. In a classical paper, Guderley and Yoshihara (1) have obtained 
(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959] 
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the solution for a wedge in a sonic stream up to the limiting characteristic 
emanating from the shoulder. Further continuation of the solution can 
then be undertaken by the method of characteristics. The given solution 
up to the limiting characteristic is complicated and some approximation 
is unavoidable in the analysis. Cole (2) has used a model in which the 
sonic line is straight and at right angles to the direction of the main 
stream. But there is no way of continuing this solution into the super- 
sonic region beyond the shoulder of the wedge. It is usually argued in 
this case that the upstream influence of the flow in this region is bound 
to be small and that the flow upstream from the shoulder may be well 
represented by a solution which ignores the influence of the region down- 
stream. More recently Helliwell and Mackie (3) discussed a solution based 
on the modification of the usual Helmholtz flow which was first suggested 
by Roshko (4) for some problems in incompressible flow. In this model 
a free streamline at sonic velocity leaves the shoulder of the wedge and 
after a finite distance becomes parallel to the direction of the main stream. 
The velocity on this streamline then decreases to the velocity at infinity. 
This solution possesses the advantage of mapping the entire physical plane 
(outside the wake) onto the hodograph plane but it cannot be said to give 
an altogether physically satisfactory description of the downstream flow. 
The object of the present paper is to formulate and solve a problem in 
which the flow closes up completely behind a symmetric body placed in 
a uniform stream with a high subsonic velocity at infinity. The body is 
to be symmetric about both axes and the Mach number is everywhere 
less than unity except on a part of the body where it attains this value. 
To enable us to formulate the problem in the hodograph plane we shall 
choose a body such that in the given flow the magnitude or the direction 
of the velocity is constant along its profile. We must bear in mind that 
the most obvious shape to choose, the double wedge or rhombus, cannot 
be considered because the sharp corners give rise to supersonic regions 
and these lead to shock waves. (An approximate treatment of flow past 
this body is given by Trilling (5).) Consequently we are led to consider the 
profile shown in Fig. 1 in which BC and EF are straight and the curvature 
of CDE is adjusted so that the velocity will be exactly sonic along it. 
The precise shape of CDE will be determined a posteriori from the solution. 
This model is suggested by the work of Riabouchinsky (6) on cavities in 
incompressible flow. 

The first step in finding the solution is the identification of the singu- 
larity in the hodograph plane. It is shown in the next section that the 
character of the singularity is essentially different from that in the solution 
discussed by Cole (2) and in the earlier paper of Helliwell and Mackie (3). 
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[t is pointed out that the particular singular solution obtained is itself a 
formal solution of the whole problem for the fictitious gas discussed by 
Germain and Liger (7). However, in order to complete the solution in 
terms of the “Tricomi gas’ it is necessary to satisfy an additional boundary 
condition on the line v = vy, the mapping of the wedge face in the hodo- 
graph plane. Inthe problems mentioned at the beginning of this paragraph 
this is achieved without difficulty by superposing an infinite system of 
images. This method does not succeed in the present case because of the 
mixed type boundary conditions on the axis v = 0. The addition of an 


s 
? 
' 
' 


—_ — 


. 1. Physical plane. 


image system, however, does lead to a solution which satisfies the boundary 
conditions on v = v, and on part of the axis. The difference between this 
and the exact solution of the boundary value problem is examined 
detail in section 3 where a complete solution is obtained as a series in powers 
of the transonic similarity parameter. It appears that the error involved 
when this difference is neglected is very small and that in most cases it 
lies well within the error inherent in the use of the transonic approxima- 
tion. Finally, in section 4, some qualitative and also some numerical 
results are given to show the type of profile for which a flow pattern of 
this nature exists. 


2. The boundary value problem 
We return now to equation (1) where uw and v are dimensionless ex- 
pressions for the velocity perturbations on a uniform flow with sonic 
velocity a*. Specifically 
(y+1)(a*—U) whee (y+)V 


u AEE 4. soa 
a* 


; = 


a* 


where U, V are the components of the actual velocity, y is the position 
coordinate measured at right angles to the uniform stream and within 
the limits of the transonic approximation is proportional to the stream 
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function ys, and y is the ratio of specific heats of the gas. The position 
coordinates x and y are related through the differential equations 


a) Ox a) 
me —_ = u y 


Ox 
év ou ~ ow ov 

It is convenient at this stage to introduce a new variable r defined by 
r = 3u'. The solution of equation (1) obtained by the separation of 


variables is then 
y= rtet\@_ (Ar), 


where ©, ,(z) is any linear combination of Bessel functions of order + 4. 
The boundary value problem appropriate to Fig. 1 can now be formu- 
lated. It is clear from symmetry that only one-quarter of the physical 


"| 





/ 


Fic. 2. Hodograph plane. 








00 


plane need be considered. Thus the upper left section of Fig. 1 is mapped 
into the semi-infinite rectangle r > 0, 0 < v < v, in the hodograph plane 
which is shown in Fig. 2. Here vy = (y+1)6 and this is the value of v on 
the straight part of the profile if 5 is the semi-angle at the nose. The 
following boundary conditions in the hodograph plane have to be satisfied: 

(i) y must have a certain singularity at O’O, the point in the hodograph 
plane which corresponds to the point at infinity in the physical plane. 


This point has coordinates (r,,0) where r, is the value of r corresponding 
to the (subsonic) flow at infinity. 


(ii) y= 0Oon OB. 
(iii) y= O0on BC. 
(iv) y=0Oon CD. 
(v) 2 = constant on DO’. 


We have used r > © to represent the stagnation condition at the front 
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of the body in accordance with the usual approximating procedure. 
Conditions (ii), (iii), and (iv) have to be satisfied since the line OBCD in 
Fig. 1 is part of the dividing streamline. (v) is essentially the condition 
which has to be satisfied if the flow is to be symmetric about a line drawn 
through the centre of the body parallel to the y-axis. It could be written 
alternatively as éy/év = 0 but it will be seen later that the present form 
is preferable because it avoids certain difficulties of convergence. 

A typical streamline is sketched in Fig. 2 and it is observed to meet the 
r-axis at right angles. If the streamline were extended across the r-axis 
it would be ‘heart-shaped’ as distinct from the double oval shape which 
is the characteristic type encountered in the singularities obtained by 
Cole, and by Helliwell and Mackie, in the problems mentioned in the 
previous section. This distinction is important. It leads to the usual 
mathematical difficulties which result from boundary conditions of mixed 
type. 

We shall now examine the singularity at O’O in more detail. It is 
apparent that y varies from 0 to 0 as r varies from 0 to r, on v = 0, that 
is along DO’. To examine the strength of the singularity we consider a 
typical incompressible problem. In this case the leading terms in the 
velocity for the flow past a body with no circulation are given by 
U-iW =U, = 5) 
(w-+-2y)* 


where U, is the velocity at infinity and C is a constant. Thus on x = 0, 


OD Air emma as U > U.. 
d (U—U)! v 1 


For the high-speed case we may assume that far away from the body the 
effects of compressibility will be negligible and we therefore look for a 
singularity of the same strength. 


In order to obtain a representation of the singularity we try a solution 
of the type 


ow 


y, =r [ f(r) Av J, (Ar) da. 


0 
The choice of J,(Ar) is made in order to satisfy the condition (iv) which 
requires y to vanish on r = 0. 


The corresponding value of x, is found 
from equations (2) to be 


x, = (§)'r8 [ f(AjeI_,(Ar) dd. 
3 
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For conditions (ii) and (v) to be satisfied it is necessary that 


[ fA)AAr) dA=0 (r >1), (3) 


0 
ri [10 ;(Ar) dA = constant (0 <r <1). (4) 


Finally, { f(A)J,(Ar) dA must become infinite like (r,—r)-*? as r>r,—0. 


It is ale verified from a suitable table of Hankel transforms (e.g. 
Erdélyi (8)), that if fd) = I (Ar,), 


then the resulting integrals satisfy the conditions (3) and (4). Accordingly 


y= rt | Ale”, (Ar, )J,(Ar) dA, 


a, = (3) i Ate” J, (Ar, )J_,(Ar) dA, 


and on v 21 


Ye = a + (0<r<n), 


= @ > > 1), 
nk 
ee ts 

In addition it may be shown that x, > 0 as r > oo which fixes the origin 
of coordinates at the stagnation point at the leading edge of the body. 

Every condition except (iii), the condition on BC, has now been satisfied. 
It is therefore necessary to add non-singular solutions to y, to satisfy the 
condition on BC but without upsetting any of the other conditions which 
have already been satisfied. It is this problem that is particularly difficult 
because of the mixed boundary conditions on the r-axis. If the boundary 
condition on v = 0 were y = 0, it would be possible to satisfy the con- 
ditions on v = 0 and v = v, simultaneously by adding image singularities 
at the points r= r,, v = 2kvy, where k = +1, +2.,..., the singularities 
being of the same strength as the original and having the same sign as k. 
When added together these would give 


while (09<r<n). 


y = 9 f ye =) Hrs) Or) dd. (5) 


For the present problem, however, this solution would fail to satisfy 
the condition on DO’. Alternatively we might attempt to add images at 
the same points but with negative signs for odd values of k and positive 
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signs for even values of k. These would give a solution 


x sinh (Vp = ”) F,(dr,d,(Ar) dd. 


1 4 
y ; cosh Av, 


0 
Now the boundary condition on DO’ is satisfied but not that on OB. 


Nevertheless, it is reasonable to suppose that the discrepancy in either 





case will be small and will decrease as vy increases compared with r, (that 
is as the transonic similarity parameter becomes small). For then the 
image singularities are further away from the r-axis and the influence 
they have in upsetting the boundary conditions on this axis is correspond- 
ingly reduced. In the following section we examine the correction terms 
which have to be added to complete the solution. 

Before proceeding to this, however, we point out that in one sense y, is 
itself a solution of the problem. This is when it is considered in conjunction 
with the transformation described briefly by Germain and Liger (7) and in 
detail by Liger (9). By means of this transformation, the hodograph equa- 
tion for a fictitious gas whose behaviour is very similar to that of a poly- 
tropic gas is transformed into Tricomi’s equation. At the same time the 
strip in which the solution of a problem of the type considered is to be 
found is mapped into an entire quarter-plane in the variables of this 
equation. It follows that if we have to solve a problem in a strip for the 
new gas, we have to solve the equivalent problem in the quarter-plane 
for Tricomi’s equation, effectively ignoring the boundary condition on 
v= Uy. Thus y, is the solution, for this new gas, of the original problem. 
By the suitable choice of parameters arising in the transformation we may 
interpret this as a flow past a thin body or alternatively we may take 6 = 47 
and consider the flow past a flat-nosed obstacle. However, the mathematical 
details of the transformation are somewhat intricate and will not be pur- 
sued further here. 

3. Investigation of the correction terms 

The solution of the boundary value problem formulated in the previous 
section will now be completed by the addition of terms to the expression 
(5) which will correct the boundary condition on v = 0 without disturbing 
those already satisfied on the other sides of the rectangle BCDB in the 


hodograph plane in Fig. 2. Accordingly the complete solution will be 
written 


y=r | {A(Ar,)*J;(Ar,)+ F(a) BB Ao—) 5 rp) dd, 
J sinh Avy 
0 x 
ax = (8)trt | {A(Ar,), (Ar) ++ FO) a) 5 oe) a, (6) 





7 sinh Avy 
0 
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where A is a constant which has been introduced as a scale factor in order 
to make the length of the straight part of the profile unity. A factor r} 
has also been introduced for convenience later. Clearly y = 0, as required 
on r = 0 and v = vp, and in addition the singularity at r = r,, v = 0 is 
of the correct strength, provided the contribution from the integral in- 
volving F(A) is non-singular at the point. This will be verified from the 
subsequent solution. 
If the conditions on v = 0 are to be satisfied, 


(3)tri [ {A (Ar,)'A(Ar,)+ F(A)}coth Avy J_,(Ar) dA = constant 
0 (0<r<r,), (7) 


ri j {A(Ar,)§Jj(Ar,)+ FA) Ar) dA = 0 (r > 1). (8) 


We recall now the relations satisfied by x, and y, on v = 0, namely, 
° 

($)ir' 
0 


A(Ar,)'j(Ar,)J_,(Ar) dA = constant (0 <r <1), (9) 





rt [ A(Ar,)4j (Ary) (Ar) dA = 0 (r > ry). (10) 
0 
If (9), (10) are now subtracted from (7), (8) respectively, a pair of dual 
integral equations is obtained for the unknown function F(A). In order 
to bring these to a form in which they can be solved, the first one is 
differentiated with respect to r. The following pair of equations is then 
obtained for F(A): 


{ AF (A)coth Avy J,(Ar) dA = — Anh (coth Avy—1)4(Ar,)J,(Ar) dA 
‘ (0<r<n), (11) 


{ F(A)A(Ar)dA=0 (r > 1). (12) 

0 
The above manipulations have been necessary in order to avoid the 
occurrence of divergent integrals. It is possible to obtain the complete 
solution by working from the beginning with a pair of dual integral 
equations derived from the boundary conditions on v = 0. The procedure 
is to regard the flow first of all as contained in a channel of width 2k and 
then to let k tend to infinity at a suitable stage in the subsequent analysis. 
The expression (5) which incorporates the correct singularity at O’O then 
appears as the first approxir ‘ion in an iterative solution of the integral 
equations obtained in this wa,.. However, this approach involves a certain 
amount of formal manipulation with divergent integrals. The reason for 
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this is that integrals representing @y/@v are divergent when v is put equal 
to zero although their limits as v > 0 exist. It follows that if the problem 
is reduced at the very beginning to that of solving a pair of dual integral 
equations obtained by specifying the boundary conditions on v = 0, the 
interpretation of certain integrals as limits when v—0 is no longer 
possible. The method which led to the equations (11) and (12) was there- 
fore adopted in order to set up the problem in a formally correct manner. 

If the substitutions r = pr,, t = Ar,, and ¢(t) = F(A) are made, equa- 
tions (11) and (12) become 


x 


ts(tjeoth("*).j (ot dt A t! othe") — WE AOEACD dt 

"7 J ry 

0 0 (0 go 8 
( d(t)J,(pt) dt 0) (p > 1). 


0 
Integral equations of this and of a more general type have been investi- 
gated by B. Noble to whom the authors are indebted for supplying them 
with the results of work as yet unpublished. By means of these results 
it may be shown that the above system can be reduced to 


[ td(t)J;(at) dt = H(x) (0<a< 1), (13) 





0 (x > 1), (14) 
where : 
H(x)+ [ AH(A)K(A, a) dA = Ga), 
K(A,a) = | ‘SxP\— "ro D) J(at) H(A) dt, (15) 
sinh(tv,/7;) 
0 
: A2?t , r ,; — r exp(—tv,/r,) 
and G j #(2— p?)-4 RZ (OS. (6) — 8 t = Gide. 
” ” i pis i ed BA (OA (08) sinh(tv,/r,) — 
0 0 


If the orders of integration in the expression for G(«) are inverted we 
find first that the integration with respect to p gives essentially Sonine’s 
first integral and, finally, that 

G(x) —AK(l,«). 


¢(t) is now found by an application of the Hankel inversion theorem to 
(13) and (14) and, if we replace the original variables, we obtain 


1 


F(A) = (Ary) [ oF (a) (Ar, x) da, (16) 


~ 


0 
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where H(«) is the solution of the integral equation 


1 
H(a)+ [ AH(A)K(A,«) dA = —AK(1, a), 
) 
and the kernel K(A, «) is given by (15). 

It is now necessary to rewrite the expression for K(A,«) as an expan- 
sion in a power series. It will be recalled that the extra correction terms 
are expected to be small when r, is small compared with v,. Accordingly 
we expand the term cosech(tv,/r,) which appears in the integrand in powers 
of exp(—tv,/r,). If the orders of summation and integration are exchanged 
it follows that 


@ 


K(A, a) = 2 [ texp(—2tr, nr, )dy(at) J, (At) dt. 
n=15 


Each term in this series may be written as an infinite sum of hyper- 


geometric functions (Erdélyi (8)) to give, after a certain amount of 
algebra, 


K(A, x) = (Aa)! s Ay (1/Vq)'V3+2™A2™, F, (—m, —§—m; ¥; a2/A*), (18) 
m=0 
(—1)"P'G+m)C(y+ 2m) 
Ww here an —_ Q2m TAT (ym! 





and ¢(z) is the Riemann zeta function. The values of a,, for m = 0, 1,..., 8 
are shown in Table | and the rapid decrease in their magnitude is apparent. 


TABLE 1 





a 


2 3 4 5 6 7 8 











0°05495 |—0°01977 | 0°006583 | —o-002085 | 0:0006371 |—0*0001896 | 000005530 























The hypergeometric functions occurring in (18) reduce to polynomials 
of degree m. If K(A,«) is expanded in powers of « we obtain 


K(A,a) = (ry/v9)™9(Aa)t 5 Byy(oe/A)2™, 
m=0 


where the 5,, are certain functions of Ar,/v, consisting of series of even 
powers starting with the term (Ar,/v9)””. 

It is now possible to obtain an iterative solution of the integral equation 
(17). A sequence H,(«) is determined by means of the relation 


1 
H, (a) = —AK(1,a)— i) AH, _,(A)K(A, a) da. 
0 
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H,(x) is taken to be 0. Then 


H,(«) = —AK(1, «) 
—A(r, Vp) N88 (b, + b, x?+-...), 
where 5,, is the value of b,, when A = 1. The next approximation is 
1 
H,(x) = —AK(1,«)+A | AK(1,A)A(A, a) da. 


0 
K(1,A) and A(A,«) are now expanded as power series and the resulting 
integration performed with respect to A. It is found that 
H,( x) —A r# (Cy +Cy x2 +Cy vft-...), 


where the c,, are certain functions of r,/v,. It is not necessary to proceed 
with further iteration as the convergence is extremely rapid and the above 
solution for H,(«) gives H(«) up to terms of order (r,/v))'®. Hence we may 
write 
A & a. oh 
H(a) = —A } ¢,, af +2m, (19) 
m=0 
Finally, we return to the original function F(A) which is given by (16). 
If the series for H(«) given by (19) is substituted in this integral it is found 
Ari 
that expressions of the form | t''®+?”"J,(t)dt have to be evaluated. On 
0 
integration by parts this becomes 


(Ar, 188+?" J, 6(Ar,) — 2m(Ar,) #2" Siai¢(Ar) 4 
}- 2m(2m— 2)(Ar,)~#+2"Jogig(Ar,)+..- « 


ence combining a 1e terms we fin 
Hence coml ll the t find 


F(A) A 5 ( ] me = Jj + m(Ary) 
— (Ar,)™ ~4 


m=1 

where the e,, are functions of 7,/v). It is convenient to define e, = 1 so 
that equation (6) can now be written 

2 = = —— 

? sinh A(v,—t 

y=A > { lye, —_ LA(Uo- ) 

= (Ar,)"-? sinh Av, 





Jy (Ar) Jy. m(Ary) dA. (20) 


+m 


) 
m = ( 0 


This expression is the complete solution of the boundary value problem. 
Table 2 gives the values of e,, corresponding to a range of values of r,/vp, 
correct to three decimal places. 


4. Analysis of the solution 

The expression (20) represents the solution for y in terms of the velocity 
variables of the flow past a thin doubly symmetric body. The velocity is 
nowhere supersonic in the flow field but sonic velocity is attained along 
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a segment of the boundary. An important physical property apparent 
from symmetry is that the total drag on the body is zero. It is now 
necessary to investigate (20) and allied equations in detail in order to find 
the geometrical ratios which determine the shape of the profiles which are 
obtained. 

First we note that for v = 0, (20) gives for y the series 
A > —1)e,, | ri(Ar,)t-" J, (Ar) Jj 


m 0 


(Ar,) da, 


+m 


and the evaluation of the integrals gives 


~* ° Dh—m,F 
Y A S (—1)*e,, — — pert 
at, (m+ $)r7 
0 (r > 1,). 
TABLE 2 


(r?—r?)"-* (O<r<n,), 
































This shows that the singularity at r = r,, v = 0 is still of the correct 
strength and that the asymptotic behaviour of y near it is not affected 
by the addition of the extra terms to the original solution (that is, the 
terms for which m > 0). 

By means of equations (2) the value of x for the complete solution can 
be obtained. It is 


=. 


—1)™e,, [ ri(Ar,)!-™ 

0 
A is now determined from the condition that the straight segment (BC 
in Fig. 1) is of unit length. We therefore let z = 1 when r = 0, v = 0. 


If r > 0 in (21) when v = v, we obtain 


cosh A(v,—v) 
sinh Av, 


Ji (Ar )I_4(Ar) ax. (21) 


~ 2% * cosech 
1= A)! > (-1)a rt [ 
(3) .\ ) ™T\(4) 1 . 

n= 0 


If the substitution t = Avy is made, the Bessel function expanded as 
a series and the orders of integration and summation of this series are 
interchanged, the above expression becomes 


Sat tea) dd. 


A6t ~ (—1)™e,, 3 (—1)"(r,/v9)?"*# fines 
n Low 22"+4n! T(m+n+) J sinht 
3 





oe 1 om 
veld) in =0 
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The integrals in this equation can be written in terms of the Riemann zeta 
function and a form suitable for computation is finally given by 





k ol « se " a mT 5 )-2n—8\r< 5\/» \$ 
bo 82S (they 5 (Tan a2 Nea A 
* ~ Ty) / ym )2ny! 1? i 111 . 
rt A '(}) —, 2 komt 22"! P(m+n+) Uo 


In Table 3 the value of Arf/v} is given for a range of values of r,/v, 
It is worth pointing out that the values obtained by using only the term 


» in the summation with respect to m contribute 99 per cent or more of 
the total. This lends considerable support to the argument put forward 
in section 2 that the correction terms would probably be of very small 
magnitude, 


TABLE 3 











542 | 09740 1°OOoI 1°036 | 1°077 








Series expansions can be found for x on the dividing streamline for all 
values of r > r,. When v = 0 in (21) 


x 


a = A(})' ¥ (—1)"e,, ( ri(Ar,)*-™ coth Avg J_y(Ar) Jj 4m (Ary) dA. (22) 


If the partial fraction expansion 


] me Av 
coth Ar, 19S o 


' \2 2 
Avy set (Av)? + (nz) 
is substituted in (22), the resulting integrals can be evaluated and we 
obtain 
oa ae, - \}—m . ( A 
nzr,\* nrr,\ ,- [nar 
. DFVWEL,.— 1,2 ™~ n 1 1 
L A2 vw Yo ; (—1)™e,, > | ) hon( 1) ("=") 
sant ae 0 | % 


(r>r,). (23) 
This gives x on the axis of symmetry as a function of r. When v = vp, 
coth Av, is replaced by cosechAv, in (22) and this has the effect of intro- 
ducing a term (—1)” in the inner series. In this way a series for x is found 
on the straight part of the profile up to the value r = r,. Bearing in mind 
the extremely rapid convergence of the summation with respect to m it 
is seen that the series provides the best way of computing the velocity on 
the dividing streamline, provided r > r, and is not too close to r,. This 
in turn leads to a knowledge of the pressure distribution on part of the 
profile. These series also serve to confirm at once the location of the 
origin of coordinates at the stagnation point for it is clear that 2 — 0 in 
(23) as r > 0. 


After the evaluation of the constant A, two other quantities have to be 
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found before the shape of the profile can be determined quantitatively. 
These are 2h, the total length of the body, and 2d, its width at the mid- 
point. This is the maximum width of the body since the free streamline 
is known to be convex to the flow. In theory a series of points should be 
plotted along the free streamline to determine its shape accurately, but 
in practice the body is thin enough for an excellent estimate to be obtained 
from a knowledge of d and h. The value of A is given by the value of x 
when v = 0 and is independent of r for r < r,. However, for facility in 
calculation, the value r = 0 is naturally selected. This will give the x 
coordinate of the point D in Fig. 1. When v = 0 and r is allowed to tend 
to zero in (21) we obtain 


r 


6+ Ss ri 
A ‘ (—1)™e,,r3-™ a m+ coth Avg Jj .,(Ar,) dA. 
Tr 1) haut 0%i+m 
(3 m=0 0 
In order to express this in a form in which h can be determined nu- 
merically coth Av, is written as the series 


cothAv, = 1+2 ¥ exp(—2kar,). 
0 - I ) 


The resulting integrals are evaluated term by term, each term containing 


TABLE 4 





























a hypergeometric function which must itself be expanded. Further re- 
arrangement of terms and summation leads eventually to the form 


Ar}(3)! S (—1)"e ee t i S (— 1" + 2n)0(3+ 2n) (2 -_ 


vy ES) n+ Bat \v) J 


—, 2™-im! a. Qhem oy (m+n+)n! 
Despite the complicated appearance of the above series, convergence is 
fairly rapid and in Table 4 the value of h is given for a range of values 
of r,/Up. 

There remains the problem of finding d. This cannot be found directly 
from the coordinates because the procedure of the transonic approxima- 
tion has identified the body contour with the line y = 0. However, it is 
not difficult to see from Fig. 1 that d is given by 


D 
= b+ | 0a, 
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where @ is the angle which the velocity vector makes with the x-axis and 
the integration is carried out along the body. An equivalent form is 

D 0 
: *vdz ‘ vx |P adv 
d o- o- “=> -- 9 


y+l c J yt! 
Cc vo 


on integrating by parts, so that finally 


d x dv. 
‘ l * 
0 
In this integral x is given as a function of v by letting r— 0 in (21). 
Hence, after performing the integration with respect to v, 


x 


A6! = rest 
A. ) ys (—1)™ 7 r m | A-% Is sm (Ary) ax. 
(y+1)0(4) — 4 
¢ o'’m=0 0 


These integrals can be evaluated without difficulty and the final form is 
d Art 3 =(*2) z = (—1)"e,, , 


r . Q 5 a | 
) V4 71 Uy a4 1.3.5...(2m-+- 1) 


m 


Table 5 gives values of d 5 for the range of values of r,/v, considered. 
































TABLE 5 
| o*4 c orf 
- ———— } 
8 j x 8 2°40 2°03 181 1°66 
PABLE 6 
| | | 
> | | 2 2 o4 e o-f 
W h rat | ” 417 0°493 =<? 602 
pe see - | 
| th rat oI 0186 0276 68 $52 
V 1)8}5 148 | 587 o7I1 Roc 0932 

















It is now possible to present the essential results collectively. A given 
value of 7,/vg determines 5d, which is the ratio of the width of the body 
at the end of the straight section to the maximum width of the body. 
We shall call this the width ratio. It also determines h, the ratio of the 
semi-length of the body to the length of a straight section. The reciprocal 
of this will be called the length ratio which is therefore the length of a 
straight section when the total length of the body is taken to be 2 units. 
In Table 6 the width and length ratios are given for different values of 
r,/U9. The corresponding value of the transonic similarity parameter 
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(1—M}) {(y+1)8}! is also given for convenience of reference. Here M, is 
the Mach number at infinity and y is taken to be 1-4. 

As might have been anticipated from the transonic similarity theory 
the width and length ratios are dependent only on the values of the 
transonic similarity parameter. However, the ratio of the width to the 
length of the body is determined only when M, and 6 are known separately. 


en eee 


Fic. 3. Typical profile for M, 0-82. 


A typical profile is shown to scale in Fig. 3. It corresponds to a (total) 
wedge angle of 10° in a flow for which the Mach number at infinity is 0-82. 
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SUMMARY 


If the cross-section of the members of a frame is varied continuously so that the 
‘strength’ of the frame follows the bending moment diagram, then, provided the 
design is carried out in a certain way, the frame will have absolute minimum material 
consumption. This theorem is proved separately here, although for plastic design 
the theorem is a special case of the general principles given by Drucker and Shield 
(1, 2). However, this special theorem applies equally well to certain completely 
elastic structures. 

The theorem is applied to determine the minimum weight of some simple con- 
tinuous beam systems, and a trial and error ‘relaxation’ technique is developed. 
These designs are compared with more practical designs in which the cross-section 
of the members varies discontinuously. 


1. Introduction 

THE calculation of absolute minimum weight of a structure is of interest 
to the designer, since he is able to estimate penalties in weight and cost 
should he depart, for practical reasons, from this minimum. It is incon- 
venient, for example, to vary the cross-section of steel members con- 
tinuously, and abrupt changes in cross-section may imply the use of only 
a few per cent more material. [Simple comparative examples are given 
below.| Reinforced concrete may be ‘shaped’ more easily than steelwork; 
here again, the designer can balance the cost of complicated formwork 
against the saving in material for the minimum weight design. 

The term ‘strength’ of a member will be left undefined for the moment. 
For the purpose of the following arguments, however, it may be of help 
to think of the strength of a member as the full plastic moment (for plastic 
design), or as the limiting elastic moment in which the yield stress is just 
reached in the outer fibres (elastic design). A more detailed discussion of 
strength, and the way in which a given strength can be achieved prac- 
tically, is given below. 

+ The results presented in this paper were obtained in the course of research sponsored 


by the Office of Naval Research under Contract Nonr-562(10) with Brown University. 
{ Fellow of Peterhouse, Lecturer in Engineering, University of Cambridge. Visiting Pro- 


fessor, Brown University, 1957-8. 


(Quart. Journ. Mech, and Applied Math., Vol. XII, Pt. 3, 1959] 
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2. Basic principles 

For a frame of given dimensions acted upon by given (fixed) loads, it 
is possible (if the frame is redundant) to construct infinitely many equi- 
librium bending moment distributions. Suppose that for one such distri- 
bution the bending moment at a generic cross-section is M,, and that the 
size of the member at that cross-section is such that the strength is also 
M,|. Although strength has not been defined, it will be assumed that if 
a frame has strength |.M,| (variable), then the total weight of the frame 
is given, to some scale, by 


W = | M,| dz, (1) 


where the integration extends over the length of all members. The impli- 
cations of equation (1) will be discussed later. 

A given function |M,| will be called a design of the frame. Of the 
infinitely many designs corresponding to the infinitely many equilibrium 
bending moment distributions, one will have minimum weight. This 
minimum weight design will be denoted |M9}. 

Consider now a deflected form of the structure defined by displacements 
y from the original position. The curvature « of any cross-section is given 
as usual by the second derivative of deflexion with respect to distance 
along the member, i.e. « = d®y/dx?. A deflected form y will be called 
compatible with a design |M,| if sgn M, = sgn« everywhere; that is, if 
sagging bending moments are denoted positive, the curvature is positive 
at a cross-section where a positive moment acts, and negative at a section 
of hogging bending moment. 

The following theorem holds: Jf, for the given loads, a design |M%\ can 
be found that is compatible with a deflected form of the structure having 
constant absolute curvature, then the design |M§\ has minimum weight. 

This theorem is a special case of the general theorem of Drucker and 
Shield (1, 2) for the continuum. However, the general theorem applies 
only to plastic structures; the following proof of the special theorem is 
valid for elastic as well as for plastic frames. 

Suppose concentrated loads P and distributed loads p are applied to 
the frame, these being in equilibrium with the set of bending moments 
Mj. Further, suppose that for this design |M}| the deflexions y of the 
frame are compatible with curvatures +x«, where « is a constant. By 
virtual work, 


W = | pydet ¥ Py = [ Mix dz. (2) 
Now, since (sgn M$) = (sgnx), the quantity W is given from (2) as 
W = f \M| |x| dar = |x| f \M*| dx = |\x|W*, 
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using equation (1). Thus the weight W* of the putative minimum weight 


design is W/ |x. 


Consider any other design ||, where, of course, M, is in equilibrium 


yy 





iT 


(a) 





(c) 








Yt: 











with the loads P and p. By virtual 

work again, 

W ( py dx+ > Py = ( M,« dx. (4) 

Thus from (3) and (4), 

Ww | 44, * da : | M,||-| dx (5) 
d K K 


since M, and « are not necessarily 


compatible. Now «/ « +1; thus 
(5) gives 
We<[|M\dex=W. (6) 


Thus W* cannot exceed the weight of 
any other design, and hence is the 
minimum weight. 

The theorem can be_ illustrated 
simply by reference to a fixed-ended 
beam, Fig. 1. Suppose the loads in 
Fig. l(a) act downwards, for simpli- 
city, so that the ‘free’ bending moment 
diagram is of the non re-entrant form 
in Fig. 1 (b). [Re-entrant bending 
moment diagrams lead to complica- 
tions; these are discussed below.| The 
infinitely many equilibrium distribu- 
tions may be constructed by super- 
imposing on this free bending moment 
diagram a straight reactant line due 
to the end fixing moments. Now the 
minimum weight theorem states that 
the beam can take up a deflected form 
with constant curvature; simple geo- 
metrical considerations show that ares 
of constant curvature can be combined 


only in the way shown in Fig. I (c). Over the central half of the beam 
there is an are of positive curvature, and the two end portions, of lengths 
1/4, have constant negative curvature. This deflexion diagram is symmetri- 
cal, and is the only one which satisfies the theorem, being completely 
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independent of the applied loads. From the theorem, the sign of the bending 
moment must be positive in the central half of the beam, and negative at 
the ends, to correspond with the signs of the curvatures. The constant 
curvature arcs must therefore join at inflexion points, and the reactant line 
is fixed as in Fig. 1 (d). The shaded diagram gives the strength of the beam 
for minimum weight; note that this strength becomes zero at the inflexion 
points, so that the final minimum weight design is statically determinate 
(consisting of a simply supported central portion of length //2 with two 
end cantilevers). 

In constructing the deflexion diagram, Fig. 1(c), it has been assumed 
that discontinuities of slope do not occur at the inflexion points. 


LemMa. The deflected form of the minimum weight design consists of arcs 
of constant curvature. These arcs must join with no discontinuity of slope. 


Suppose there is a deflected form of the minimum weight beam for which 
discontinuities in slope @ occur at the inflexion points. Consider a design 
M%* identical with |Mj| (the minimum weight design), except that in 
the neighbourhood of the inflexion points Mj* + 0. Thus the design 
M}* is too ‘strong’, and 
wees W*. (7) 
Now, as before, __ 
W = |\x|W* = |\n|W**+ > |MG*.9), (8) 
the last term representing the work done at the discontinuities. From (8), 
W** < W*, which contradicts (7). Thus the lemma is established. 


3. Simple designs 


It was seen that the minimum weight design of a fixed-ended beam is 
independent of the loading system, provided the loads all act vertically 
downwards. The same is true of a propped cantilever; geometrical con- 
siderations show that the point of inflexion must always occur at a distance 
//,2 from the prop (Fig. 2). This independence of the minimum weight 
design on the loading system holds for only the simplest problems. 

Consider, for example, the two span beam with one redundancy sketched 
in Fig. 3(a). The bending moment diagram will be of the general form 
shown in Fig. 3(6), where the single redundancy has been taken as the 
moment at the internal support. As sketched, it will be seen that there 
are two inflexion points distant a, and a, from the internal support. 
|Most loading conditions, but not all, will imply two inflexion points for 
the minimum weight design. Attention will be confined for the moment 
to the case illustrated.] From the theorem, the deflected form of the 
minimum weight design is as sketched in Fig. 3(c), with hinges at the two 
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sections defined by a, and a,. [Note that this structure now has one degree 
of freedom, compared with the just statically determinate structure of 
Fig. 1(c).] Now a, and a, are not uniquely determined from purely geo- 
metrical considerations of the deflected form. General expressions are 
given below; it may be noted here that the constant curvature condition 


gives only one relationship between a, and a,. Some other equation must 








A B A 
(a) 








= +4 (c) 


1. 3 Fia. 3 


| ee 


therefore be found to solve the problem; this equation is furnished by 
inspection of Fig. 3(b). It will be seen that a, and a, are related by the 
actual shape of the free bending moment diagram. This implies that the 
loading conditions must be used to find the minimum weight design, and 
this is no longer independent of the loads as it was for the single span 
fixed ended beam. 

[t will now be demonstrated (but not proved rigorously) that it is 
always possible to construct a proper minimum weight deflexion curve, 
i.e. one consisting of ares of constant curvature. Consider a portion of 
a continuous beam between two supports, and suppose that the free 
bending moment diagram is drawn and a reactant line superimposed. 
Now an are of constant curvature is specified by two arbitrary constants; 
if there is no inflexion point in the bending moment diagram, the known 
deflexions at the supports (zero for simplicity) will enable the two constants 
to be calculated. If the free bending moment diagram is re-entrant, so that 
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there may exist inflexion points, a similar argument may be used. Between 
one end of the beam and the first inflexion point, an are of constant 
curvature may be specified in terms of one unknown constant (the other 
constant being determined by the condition 6 = 0 at the end support). 
By continuity of slope and deflexion at the first inflexion point, the arc 
of reverse curvature between the first and second inflexion points is then 
uniquely determined in terms of the one unknown constant. Thus each 
are may be constructed until the other end of the beam is reached; at the 
other end, the condition 6 = 0 enables the unknown constant to be 
calculated. 

Suppose now one extra internal support is provided for the span con- 
sidered. One extra condition must be satisfied at this support (6 = 0, say). 
But the extra support has introduced one extra redundancy into the 
problem; this means that the reactant line can be adjusted with one degree 
of freedom in order that the extra deflexion condition is satisfied. Similarly, 
a fixed end to a beam introduces one extra geometrical condition (@ = 0), 
but one redundancy is also introduced. 

For every example of a framed structure, the degrees of freedom of the 
reactant line (i.e. the degree of redundancy) will always enable a minimum 
weight deflexion curve to be constructed. 


4. Slope-inflexion equations 
A relaxation method is developed in the next section which ‘balances’ 
angles of rotation at we For this purpose, certain elementary solu- 
ot 


Ki. tal gt igh Nhe 
-_ a ponci—ted aa 


Fic, 4 


tions are required, and the most common cases are illustrated in Fig. 4. 
In Fig. 4 (a) a single span is bent into constant curvature with no inflexion 
points; one inflexion point occurs in Fig. 4 (b), and two in Fig. 4(c). Taking 
the constant curvature « equal to unity, the following expressions may be 
obtained for the angles of rotation at the supports: 


No inflexion points: 0, = 4 
(9a) 
,=— 


: 
_ 
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; : ; l l—a)? 
One inflexion point: 6, =>5 a — 
7 (9b) 
l @ 
On — - ] 
; l l—a,)? 3 
Two inflexion points: 0, = = a 7 
“ (9c) 
t (t—a, a? 
‘a _ 4 _ ° 
: 2 l l 


[t will be seen that if the positions of the inflexion points are known, the 
angles at each end of a span may be calculated immediately from equa- 
tions (9). Since the positions of the inflexion points depend, in general, 
on the bending moment diagram, no general solutions are possible except 
for the trivial cases of single-span beams. [For example, a single-span 
fixed-ended beam has 6, = 6, = 0; equations (9c), assuming two in- 
flexion points, give a, = a, = 1/4.] 


5. Numerical example 


The problem shown in Fig. 5 will be solved numerically to illustrate the 
relaxation technique. A continuous beam has supports at intervals of 
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8 units, each span carrying a central concentrated load of 50 units; the 
minimum weight design is required. The problem is started by construct- 
ing the free bending moment diagram, and superimposing a guessed 
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reactant line, as in Fig. 6. For this problem, the height of the free bend- 
ing moment diagrams is 100 units, and the values of the reactant moments 
at the internal supports have each been taken as 50. These starting 
moments are shown in Table 1. From Fig. 6 the values of a, and a, (see 
Fig. 4(c)) may be read off, and the rotations 6 computed for each end of 
the separate spans by means of equations (9c). The rotations agree except 
at support B, where there is a discontinuity in slope of 1-12. The value 
of M, is adjusted to make the slope continuous at B; it will be seen that 
increasing /, from 50 to 64 causes the discontinuity to decrease from 
1-12 to —0-04. A discontinuity is introduced by this process at support C; 
moment .V/,, must be adjusted, and the work proceeds until all discontinui- 
ties are small. The final reactant line for minimum weight is shown in 
Fig. 6. 

This method can of course be used for any system of loads; the free 
bending moment diagram is drawn, and the values a, and a, defining the 
inflexion points read off for the current position of the reactant line. 


6. The ‘strength’ of frames 


The minimum weight designs proposed above can be achieved with 
either entirely elastic or completely plastic frames; in both cases it is 
assumed that a suitable load factor is incorporated in the loads. 

Examining first an entirely elastic structure, it will be seen that the 
conditions of the basic theorem will be satisfied if (1) the yield stress is 
just reached in the extreme fibres at all cross-sections, (2) the members 
have constant depth, and (3) the width of the (rectangular) cross-section 
is adjusted so that the moment of resistance at each section equals the 
design bending moment. Condition (1) is one of the strength criteria; all 
sections are fully stressed. Since all cross-sections have equal maximum 
stress, the requirement of condition (2), that all members have the same 
depth, ensures that constant curvatures are developed. Condition (3) 
completes the strength criterion and ensures that the weight of the mem- 
bers is a linear function of the moment of resistance, as implied by 
equation (1). 

These three conditions could be almost exactly fulfilled in practice with 
materials like wood, or like pre-stressed concrete; note that condition (3) 
implies that the width of the members is zero at the inflexion points, so 
that some practical adjustment must be made at these sections. The three 
conditions could be satisfied reasonably closely with an I section steel 
beam, the flanges being curtailed continuously or in discrete steps. The 
linear weight criterion of equation (1) would not be satisfied so closely 
with this type of design. 








YMAN 


. 
4 


oe 
_ 
— 


J. 


£0-0 


cO-0 








0g 


£0-0 


10°% 


tn 


CO-O 


OL-0 


Ty ‘ty ysnlpy 

















0-0 


09-0 





on ynlpy 




















9-0 | ¢ 





9-0 





# yy qyenfpy 


qIBIS 





| AIAV, 





ABSOLUTE MINIMUM WEIGHT DESIGN OF FRAMED STRUCTURES 323 


Suppose now that the material used exhibits an elastic-ideal plastic 
stress-strain curve, a condition closely approximated by structural mild 
steel. If the full plastic moment of each cross- 





section defines the ‘strength’ of that cross-sec- PY 


tion, then all sections of the structure must be Te 











fully plastic in order that the minimum weight 
theorem be satisfied. Since curvatures can 
have any values under these conditions, the 
curvature condition of the theorem will be 








satisfied automatically. Thus the only prac- 
tical problem with a fully plastic design in 
steel isthe actual variation of the cross-section 
to follow the minimum weight bending mo- 
ment diagram; the theoretical problem is to 
make this variation in such a way that the 
linear weight condition is satisfied. Again, 














these conditions can be satisfied reasonably 
closely if I beams with curtailed flanges are 
used, 


7. Practical designs 

In practice it will not be possible to vary 
the cross-section of a member continuously, 
or even discontinuously at more than a few 
sections. It is of interest, therefore, to 
examine the consequences of using prismatic 
members over some, at least, of the length 
of the members, and, in particular, to esti- 
mate the weight penalty involved in such 
modified designs. As an illustration, only the 
example of the fixed-ended beam will be 
discussed. 

Fig. 7(a) shows the collapse bending 
moment diagram if a prismatic beam is used 
to carry the uniformly distributed load p. 
M, is constant at the value pl?/16, and the 
weight, from equation (1), is 














m]3 3 
WP oP (10a) 


Fig. 7(b) shows the bending moment distribution for minimum weight 
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design, the inflexion points occurring at quarter span; , is given by 


pa SM Rha 


32 2° 2] 
leading to the minimum weight 
, pl 
W = = (10 b) 


Thus the prismatic beam uses just double the material of the minimum 
weight beam. 

A more practical design is shown in Fig. 7(c). A prismatic central 
portion of the beam connects two linearly tapering end portions. Thus, 
if the moment M in Fig. 7(c) has the value 4pl? (corresponding to the 
minimum weight bending moment distribution), the prismatic central 


portion has M, = ,,pl*, and the corresponding weight is 
3 
W = 1-29 PE (10c) 
32 


The central prismatic portion has length //v2. In fact, if M is chosen as 
5pl?/48, a relative minimum weight design is obtained for the type of 
beam sketched in Fig. 7 (c): 
, 3 , 
W = 1-232, (10 ¢) bis 
32 

Finally, Horne (3) has given the weight of a beam of the shape sketched 
in Fig. 7 (d); the relative minimum for this shape is 


a pl 
W 1-567_. (10d) 
32 
These figures give an indication of the economies which may be achieved 
if non-prismatic members are used. 
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THE GENERATION OF TORSIONAL STRESS WAVES 
IN A CIRCULAR CYLINDER 


By R. P. N. JONES+ (University of Sheffield) 
[Received 11 June 1958] 


SUMMARY 


A general solution in series form is obtained for problems in which torsional waves 
are generated in a bar of circular section as a result of given initial conditions, or 
given external stresses. The solution is in effect an extension of the Pochhammer-— 
Chree theory of harmonic wave propagation, and it satisfies exactly the elasticity 
equations and the boundary conditions. A particular problem is considered in which 
there is a suddenly applied torque in the form of a distributed shear force acting 
circumferentially on the surface of the bar. In this problem the displacements at 
the surface are given satisfactorily by the elementary theory at large distances 
along the bar, except near the wave front, where the convergence of the series is 
slow. The derived series for stress and velocity at the surface of the bar do not 
converge satisfactorily. An alternative solution in terms of waves of distortion in 
an infinite medium is discussed. 


1. Introduction 


Ix the elementary theory of torsional vibration of a uniform bar of circular 
section, it is assumed that each cross-section undergoes a pure rotation 
¢, about the axis of the bar. This assumption leads to a simple wave 
equation, viz. ed, 2p 

ar as One (1) 

ot? 2? 
in which the axis Oz lies along the axis of the bar and c? = y/p, where p 
is the density and yw the modulus of rigidity. It may be shown that this 
solution satisfies exactly the equations of motion of an elastic solid and 
the boundary conditions on the free surface of the bar. However, in the 
Pochhammer-—Chree theory of harmonic wave propagation (see Kolsky (1) 
for references) it is shown that additional types of torsional wave motion 
are possible, in which the displacements over the cross-section are of more 
complicated form, having one or more nodal circles. To determine the 
importance of these higher modes of wave-motion, and the limitations 
of the elementary theory, it is necessary to consider the generation of 
torsional disturbances by a given initial motion, or given external stresses. 
A general method of obtaining exact solutions to such problems is described 
in this paper. 

+ Now at the Engineering Department, Cambridge. 
(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959) 
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2. Solution of equation of motion 
It is convenient to use cylindrical coordinates (r, 6,2) with the z-axis 
lying along the axis of the bar. If the corresponding displacement com- 
ponents are denoted by (w,v,w), then, in the case of torsional motion, 
u and w vanish everywhere, and v is independent of 6. With these 
simplifications, the equations of motion of an isotropic elastic medium (2) 
reduce to the single equation 
ev “(= es a (2) 
ot — oe Fer 
All stress-components are zero except o9,, 0,9 which are given by the 
equations 


C9: pL ‘ (3) 


0,9 pr : (“) (4) 
cr\r 


If the surface of the cylinder is free from stress, the boundary condition 


to be satisfied is that p 
[o,9], - 0, (5) 


where a is the radius of the cross-section. 

It has been stated above that (2) and (5) are satisfied by the elementary 
solution ; 
v = rd,(z,t) = % (6) 
where ¢, satisfies (1). It may be verified that another solution of (2) 
satisfying the boundary condition (5) is 


v ad, ( K,, r)d,,(z, t) v (7) 


n 
where » denotes a positive integer, J, denotes the Bessel function of the 
first kind, ¢,, satisfies the equation 
c 2h od 


ay Yn 2Kr2 
m as € Cc ke no (8) 
ot- C2" 


and Ka is the nth non-zero root of 


= F 1(Ka)| «%, (9) 

that is, J,(Ka) = 0. (10) 

The first nine roots of (10) are given by Jahnke and Emde (3). In particular 
K,a= 5-136, K,a=8-417, K,a= 11-620, K,a~ (n+3)n 


when v is large. 
The Pochhammer-Chree theory deals with the propagation of harmonic 
waves. Solutions of this type may be obtained by putting 


dn, = explik(z—Vt)], (11) 
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where k is the wave number and V the phase velocity. The case n = 0 
corresponds to wave propagation in the fundamental mode, which is 
dispersionless (_V = c), whilst the solutions for n = 1, 2,... correspond to 
the higher modes of torsional wave propagation given by the Pochhammer-— 
Chree theory, in which 


(12 


It will be assumed that the general solution of (2) satisfying the boundary 
condition (5) may be expressed as an expansion in terms of the charac- 
teristic solutions (6) and (7), i.e. 


«2 


= z Un _ rpy+a > $,, J,(K,, r), (13) 


n=0 =1 


where ¢, and ¢,, satisfy (1) and (8). This general solution may be regarded 
as an extension of the Pochhammer-—Chree theory of torsional wave 
propagation. The first term of the series is of course the solution given 
by the elementary theory. 


3. Initial value problem for bar of infinite length 


Consider a bar of infinite length in which the initial displacement and 
velocity are given by the equations 


v = f,(r,z) 


év . (14) 
= = f,(r,2) 


Assuming that (13) may be differentiated term by term, the initial condi- 
tions (14) give 


filrs2) = rdol2,0) +0 ¥ A(Kyr)b,(2,0), (15) 


falrs2) = rbol2,0)-+4 Y J(Ky Mibn(2), (16) 


where dots denote differentiation with respect to t. It may be shown that 
the characteristic functions J,(K, 17) satisfy the following orthogonal 
properties : 


[ AUK, rr? dr = 0, (17) 
0 


[ A(K, r)(K_r)rdr=0 (n #m). (18) 
6 
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By use of these properties equations (15) and (16) may be solved for 
¢,(z, 0), ete. Thus 


2,0) = 5 | f(r, z)r? dr, (19) 


0 
rf. 
¢,,(z,0) = Pa | flr, 2)J,(K,, r)r dr, 


a 


where rT. | [A(K,, r)Prdr = 4f[aJ,(K,,a)]?. (21) 


. 


0 
Similar expressions are obtained for ¢,(z, 0), ¢,,(z,0). The solution is com- 


pleted by solving (1) and (8) with the given initial-value functions ¢,(z, 0), 
ete. The solutions are 


z+cet 
* 


l 
ct, 0)+44,(z—et, 0) 4 : | o(€, 0) dé, 
oC 


: , os —— 
J\(K,, 0)d,(€, 0) dé +-— | J,(K,, 0)¢,,(€,9) dé, (23 
a | 
z—ct 
where o? = (ct)?—(z—€)?. 


4. Bars of semi-infinite or finite length 


When the bar does not extend from z = —o toz = +-0 it is necessary 
to satisfy additional boundary conditions on the ends of the bar. In the 
simple cases where the end of the bar is either ‘free’ (a9, 0), or ‘fixed’ 
(vy = 0), solutions satisfying the boundary conditions may be obtained 
from the infinite bar solution by the method of images. If, in the infinite 
bar problem, the initial displacement and velocity are even functions of z, 
then considerations of symmetry show that there is no stress at the section 
z = 0 (unless an external torque is applied at that section) and the condi- 
tions at z = 0 are those of a free end. Similarly, if the initial displacement 
and velocity are odd functions of z, the conditions at z = 0 are those of 
a fixed end. A similar method may be used to discuss the reflections of 
waves at a fixed or free end, and in either case it may be shown that the 
reflected wave is of the same form as the incident wave. At a free end 
the displacements of the reflected and incident waves are of the same 
sign; at a fixed end they are of opposite sign. 

An important problem is that in which stress waves are generated in 
a semi-infinite bar (z > 0) by the application of shear stress og, on the 
end z = 0. Let the given applied stress be 


[92]. o = f,(r,t). (24) 
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Assuming ee of (13) we have 
hal r,t) = rdo(0,t)+a > J(K,, r)¢,(0, t), (25) 


where primes poten differentiation with respect to z. By use of (17) and 
(18) we obtain 


4 a 
$,(0, t =f tyr? d 26 
0 ) pat J 3(", rs ar, (26) 


¢,,(0,t) = —— ar | t)J,(K,, r)r dr. 


If ¢, and ¢,, are initially zero every ae the solutions of (1) and (8) are, 


for ct < , 
aiid do — dy, =0 
t—zic 


$o(z,t) = —e | (0,7) dr, (28) 
0 


and for ct 


t—z/c 


bn(2,t) = —c | I(K,8)6,(0,7) dr, (29) 
0 


where s? = c?(t—r)?—2?. (30) 
The solution given by (28) and (29) exhibits a wave front at z = ct. 
By differentiating these equations and proceeding to the limit ct = z+-0, 


we obtain ‘i 
#,(2,2) = #4(0,0), (31) 


4,(2.5) = #4(0.0) = —e$5(0,0), (32) 


and therefore ¢’ and ¢ are constant at the wave front. ‘It follows that the 
distributions of stress and velocity across the section of the bar at the 
wave front are constant and equal to those obtaining at z = 0, t = 0, i.e, 


(opclena = Salt, 0), (33) 
[6a = ——— halt, 0). (34) 
v(HP) 


5. Bar with external torque applied at surface 
The general solution given above will now be applied to a problem 
where stress waves are generated by the sudden application of an external 
torque at the surface of the bar. Consider first a semi-infinite bar, z > 0, 
and let there be a shear stress over the section z = 0 given by the equations 

=0, 0<r<a—h 
ae = ArH(t), a—h<rc<a , 
Zz 


(35) 
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where H(t) is the unit step-function and A is a constant, given by the 
equation oT 
Pf ae ees. Sas (36) 
nlat*—(a ~h)*| 


where 7), is the resultant torque of the stresses (35). From (26) and (27) we 





obtain —_ 
4,(0,1) = —=—8, (37) 
d a7pas 
,,(0, t) = _Ala— h)PJ| K,(a—h)] (38) 


pak, P, 


By proceeding to the limit h = 0 we obtain the case where the external 
torque is applied in the form of a line shear force, of magnitude 7),/2z7a? 
per unit length, acting around the circumference r= a at z= 0. The 
solution of this problem will also apply (for z > 0) to an infinite bar with 
an external torque 27), H(t) applied at r =a, z = 0. 
In the limit / 0, (38) becomes 
- A ” 
®, (YU, t) ~ 7uad,(K,, a)” (39) 


and from (28), (29), we have, for ct > z, 





9p 
$o(2,t) = = (ct—z), (40) 
a- 
t—z/c 
Be : . 
$,(2,¢) = ——— Jo( K,,8) dr, (41) 
a*J,(K 
0 
1 
where B= —° (42) 
TA" |b 


and s is defined by (30). 
The displacement at any point on the surface of the bar is given by the 


equation : 
[ena = > [nha (43) 
n 0 
in which, for ct > z, 
: 2B 
Vol, 1 2 (ct —z), (44) 
a 
; t—2/c 
, sc 6 ; " 
[vn], —— — | J(K,, 8) dr. (45) 
0 


~  2a0)(n + })-#sin(7B) (46) 


me 


where x? = (ct)? —2?, (47) 
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and B = (n+ ) _— (45) 

é ( = (n-+ — <= =» < 
a 4 

The series (43) is absolutely convergent, although at small or moderate 


Vo 








oO 


Fic. 1. Displacement components at surface of bar for z = 4a. 





v; x10 


v, x10 
-0+1 = 


Fic. 2. Displacement components at surface of bar for z = 40a. 


values of z/a, and in the region of the principal wave front, the convergence 
is too slow for satisfactory computation. However, for large values of z/a 
the terms ?,, v,,... are small, and the displacement is given fairly accurately 
by v, only, except at the wave front. This is illustrated in Figs. 1 and 2, 
which show the displacement components %, v,, v2, at two positions on 
the surface of the bar, for z = 4a and z = 40a. The components v, and 
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v, are in each case plotted to a scale ten times as large as that of v). By 
differentiation of (44) and (45) we have 


ol (49) 
a 

| Vn |, a = = J,( K,, x). (5 ) 
a 


For large values of n 


[¢, <a ~ mal (a) +-3)- cos(7B) (51) 
7 


ax 


oa 

and therefore the series 3 5 oe (52) 
r= 

does not converge absolutely. If « + 2am, where m is any positive in- 

teger, the terms may be bracketed into groups having alternately positive 

and negative signs, and the series is conditionally convergent. If « = 2am, 

the series diverges, since all terms are of the same sign for n +0. The 


series also diverges at ct = z+-0, since 
, ' Be ‘ 
lim [¢,|,-. =—- (53) 
ct—>z+0 a 
A similar discussion applies to the series 
2 
> [vn], a’ (54) 
n=0 
The divergence of (52) and (54) at ct = z+-0 may be inferred from equations 


(33) and (34). The shear stresses og, at z = ct are the same as those applied 
at z = 0, t= 0, and consist of a line shear force acting around r = a. 
The shear stress and velocity at r = a, z = ct, are therefore infinite. 

To avoid the divergence of (52) and (54) it is necessary to allow A in 
(38) to remain finite. Alternatively, in the case of the infinite bar, the line 
shear force at r = a, z = 0, may be replaced by a shear stress 


ae te 2 
Cn-g = T, aT “€ 


distributed over a surface element r = a, |z| < ¢, in which case the solu- 
tion is obtained by integration with respect to z. Either of these modifica- 
tions is sufficient to ensure the absolute convergence of the two series. 


6. Wave solutions 


An alternative to the Pochhammer-—Chree theory for uniform bars is 
to consider the propagation of waves of dilatation and distortion, and their 
reflections at the surface of the bar. This method has been discussed 
qualitatively by Kolsky (4) for the problem of longitudinal wave-propaga- 
tion. A similar discussion may be applied to the problem considered in 
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the preceding section. In the case of torsional motion, only waves of 
distortion are present, and, because of axial symmetry, it is sufficient to 
consider wave-propagation in the plane @ = 0 only. 

In Fig. 3, P is any point on the surface of the bar (in the plane @ = 0), 
and O,, O, are points lying on the circle of application of the external 
torque r = a, z = 0. Two waves of distortion will arrive at P via the 
direct paths O, P, O, P, the displacements associated with the two waves 
being of opposite sign, since the shear forces at O, and O, are in opposite 
directions. Subsequently, reflected waves will arrive at P via paths such 








0, A; P 
“a. A Yi 
~ #% va 
\ \ 4 Tf 
{ 4 // 
ly Vv / 
/ V~ \ / f 4, 








kL Vv \ r 
0, A, B, 
Fic. 3 Fic. 4 


as O, A, P, O, A, B, P, and so on. At each point of reflection the displace- 
ments Of the incident wave are parallel to the free surface, and the reflected 
wave consists of a single wave of distortion, with angle of reflection equal 
to the angle of incidence. The arrival times of the direct and reflected 
waves are given by the equation 

a = 2ma, (55) 
where a is defined by (47), and m = 0 ora positive integer. The divergence 
of the series (52) and (54) is thus associated with the arrival of direct and 
reflected wave-fronts at the point P. 

An approximate expression for the displacements near the surface of 
the bar at the first wave front may be obtained by neglecting the curvature 
of the surface, and (in the case of the infinite bar) considering the problem 
of a semi-infinite elastic solid. Choosing new cylindrical coordinates 
(r,,9,,2,) as in Fig. 4, let the solid be bounded by the plane 6, = +4. 
A line shear force, of magnitude 7, H(t)/za? per unit length, acts along 
the line O,2z,. If the displacement components are denoted by (w,, v;, w,), 
it will be assumed that u, and v, vanish and that w, is independent of @, 
and z,; w, is then given by the single equation of motion 
oem ae 


= om 56 
et? or? 1, Or, (56) 
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All stress components vanish except o,, which is given by 


Cw 


7) 


ct 
9.9 : 
7a" P 


Ww; 


cosh- 


(ct > r,) 


ry 


satisfies (56) and also the loading condition 


T, 
. . 0 <4 
lim [—7r,o,,] = —* H(t). (59) 
r—0 7a“ 

Transforming to the original coordinates and notation, we have an approxi- 
mate expression for the displacements at the surface of the bar in the 
neighbourhood of the first wave front, viz. 


ct’ 


[v], = = cosh-*( (ct > 2). (60) 


7 ~ 
\ 


It is to be expected that this approximate wave-solution will be satis- 
factory only when the disturbance is confined to the outer layers of the 
bar, i.e. when a is small compared with a. Equation (60) therefore applies 
only in the region of the first wave front, for small or moderate values of 
z/a. In these circumstances the convergence of the series solution (43) is 
poor, and a numerical comparison of the two solutions is difficult. The 
solutions are plotted in Fig. 5 for z = 4a and values of (ct—z)/a up to 
4x 10-8. At the maximum value of the abscissa the disturbance at z = 4a 
will have spread to a radial depth of 0-18a approximately. In the case of 
the series solution three curves are given, representing the sums to one, 
four, and ten terms. In view of the approximations involved, agreement 
between the two solutions is reasonable. 


7. Conclusions 

The generalized Pochhammer-Chree theory of torsional waves given in 
this paper shows that the importance of the higher modes of wave propaga- 
tion depends greatly on the manner in which the waves are generated. 
Only in the particular case where the initial conditions or the applied 
stresses themselves conform to the elementary theory will the higher 
modes be absent. 

When torsional waves are generated by shear stresses at the end of a 
bar, the higher modes are important at the wave front, as well as at the 
end where the stresses are applied, and the elementary theory is unsatis- 
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factory in these regions. The example considered in the paper is an 
extreme case in which there are suddenly applied stresses in the form of 
line shear forces around the circumference at the end of the bar. In this 
example the displacements are given satisfactorily by the elementary 
theory at large distances along the bar, except near the wave front, but 
neither the elementary theory nor the generalized Pochhammer—Chree 
theory gives a satisfactory solution for the stresses or velocity, since the 





i —; £3 «10° 


2 3 


Fic. 5. Displacement at surface of bar for z = 4a. W = approximate 
wave solution (60). S,, = sum to m terms of series golution (43). 


series fail to converge at certain points. In this case an alternative solution 
in terms of waves of distortion is desirable. However, it is probable that 
the generalized Pochhammer-—Chree theory will be adequate for the calcu- 
lation of stress and velocity in less severe cases, where the external stresses 
are distributed over a finite area, or are applied gradually. 

There is little experimental evidence available on torsional waves, apart 
from that obtained by Davies and Owen (5), (6), with a torsional pressure- 
bar. In these experiments, torsional waves were generated in a long bar 
by the application of a torsional impact at the mid-point, and the waves 
were observed by measuring surface displacement at a considerable 
distance along the bar. These measurements showed no evidence of the 
excitation of the higher modes, a result which is in agreement with the 
conclusions of the present paper. 
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A NOTE ON THE APPROXIMATE CALCULATION OF 
THE TEMPERATURE DISTRIBUTION IN AN 
INCOMPRESSIBLE LAMINAR BOUNDARY LAYER 
OVER A HEATED PLANE SURFACE 


By D. E. BOURNE (University of Sheffield) 


[Received 24 July 1958] 


SUMMARY 
The analysis of an earlier paper (1) is completed by evaluating in terms of 
Whittaker functions approximate expressions, in the form of definite integrals, for 
the temperature distribution in a laminar boundary layer over a heated plane 
surface. 


1. Introduction 


Ir has been demonstrated by Davies and Bourne (1) that the temperature 
distribution in an incompressible laminar boundary layer over a heated 
semi-infinite plane surface can conveniently be calculated by an approxi- 
mate method based on a power law representation of the exact velocity 
profile. Consider the particular case, when the mainstream velocity U, 
is constant and the difference between the mainstream temperature 7, 
and surface temperature 7, is given by the series 
+ 
T,—T, = > 4,2", (1) 
n=0 

where x denotes distance from the leading edge of the surface and the a,, 
are constants. Then the temperature 7' in the boundary layer is given by 


T—T, = ¥ a,2"X,(é), (2) 
n=0 
where the variable € is defined by 


cin (=). (3) 


2\ va 
y denotes the normal distance from the surface, v is the kinematic coeffi- 
cient of viscosity, and the functions X,(€) are expressed as certain in- 
tegrals. In the case n = 0 the integral was shown (1) to be expressible 
in terms of a Pearson J-function (2), but for n > 0 no reduction was 
found possible. 

It is the purpose of this note to complete the analysis by showing that 
in general the integrals are expressible in terms of Whittaker functions 
(3, p. 304). 


(Quart. Journ. Mech. and Applied Math. Vol. XII, Pt. 3, 1959) 
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2. Reduction of the integrals 

Denoting the velocity component, in the boundary layer, parallel to the 
surface by U, the velocity profile is represented (1) approximately by the 
formula Viv bebe, (4) 
where 6, 8 are constants. The functions X,(&) are then given (1, equation 
47)\+ bv : . 
(4¢))7 by X,,(€) 7m tu(sin wr)(1+2n)B,{u, 1+2nphL,, (5) 
where B, is a beta function, 


m (1+-8)/(2+-B), 


| — 2b, o€2+8 
and 4 (1—w)-hyei+2n) lexp! _— 26, of*** | 


—y * 
; \(1+8)(2+8)(1—u)} 


0 


o being the Prandtl number. J, is the integral to be evaluated. 


2b, o€2+8 
(1+8)(2+-8) 
and changing the independent variable to v = nu/(1—u), the integral 
becomes 


We write (8) 


x 
. 


y\ —1-2un 

2n) pul +2n) 1 1 ‘ e-* dv. (9) 
\ ” 

Setting 


p(1+-2n) -— $41; -1—2un = k—}+H, 


l= }(u—1), (10) 


- »\ k-}H 
we have é' ety er | y-k-tHi 1] 4-— e-" dv, (11) 
: 
0 


and using the definition of the Whittaker function W,,(y) given by 
Whittaker and Watson (3, p. 340) this may be written 


I, = e784 +2—)-kD (4 —k+1)W, (9). (12) 
Substituting k, / from (10), we find from (5) and (11) that 
X (9) = w-y(sin p)(1+-2n)B, fu, 1+-2np} > 
<e-9y- eT (un + 2un)W_pion+t,au—v(%)- (13) 


The result may be simplified by expanding the beta function in terms 
of gamma functions and using the relation 7—(sinu7)I'(u) = {f(1—p)}" 


+ An error in that result is corrected here. 





A NOTE ON TEMPERATURE DISTRIBUTION 


(3, p. 239); after substituting from (6) and (8) we find 


X,,(€) = {P| . r+ een 2b,0 “1 +Pvaa+h) 


. é v-Pexpl 





“er x 


J) (© (2+8) J\U+8)(2+8)) 
—b, of? |, (4n+1)(1+8) 1 2b, ofth . 
(1+B)(2+B)) ~ 2a+B* 32+B|(1+8)(2+8) 
This is the required result. 
It is noted that the result, mentioned earlier, that the function X,(€) is 
expressible in terms of a (tabulated) Pearson J-function, follows easily 
from (14) (cf. 3, p. 341, noting that W,,,(z) = W,_,(2)). 





(14) 
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DIFFUSION OF IONS IN A STATIC F, REGION 


By J. E. C. GLIDDON 
(Queen Mary College, Mile End Road, London, E. 1) 


[Received 20 November 1958] 


SUMMARY 


The vertical diffusion of ions under the action of gravity and of a rate of electron 


} 


loss which decreases exponentially with height is shown to correspond to heat 
conduction in a rod with a heat-source distribution and radiation from the lateral 
surface. Electron density is determined as a time-periodic function of the height 
and is expressed in the form of an infinite series. 


1. Notation 
height above a reference level. 
time in seconds measured from sunrise. 
time in radians measured from sunrise. 
(= 1-37 « 104) number of seconds in a radian. 
electron density. 
number-density of neutral molecules. 
rate of ion-production. 
temperature of the atmosphere. 
coefficient of diffusion. 
scale height. 
height above reference level of maximum ion-production. 
K K(z) coefficient of attachment. 


2. Introduction 

It is now generally believed that the rate of electron loss in the F, region 
of the ionosphere is of the attachment type (that is, proportional to the 
electron density) and is also controlled by the downward diffusion of ions. 


In their discussion of this problem, Ferraro and Ozdogan (1) assumed 


that the loss coefficient, A, of the electrons was constant. This assumption 
was made to simplify the analysis and also because at that time it appeared 
that the molecular density of the air at the level of the F, region was such 
as to make diffusion the main controlling factor. However, it has been 
pointed out by Ratcliffe et al. (2) that, as was shown by Bates and Massey 
(3), under certain circumstances the loss coefficient K is not constant but 
proportional to the local molecular density of oxygen (O,) so that it 
decreases exponentially with height. The object of the present paper is 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959) 
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to give an analytical solution of the diffusion equation for the case when 
K is of the form suggested by Ratcliffe and his co-workers, namely, 
K = 10-texp( sec, (1) 
50 
The application to ionospheric problems will be reserved for a later 
paper. In particular, the calculations of Ferraro and Ozdogan will need 
some revision in the light of new satellite data which indicate that the 
density in the F, region is about 10 times greater than the value which 
had been derived from an analysis of rocket data. The new estimate of 
the molecular density is, moreover, in accord with that which had pre- 
viously been inferred on theoretical grounds (4, 5). 


3. Statement of the problem 
The equation governing diffusion of ions in the ionosphere is (1) 
aN JAN 3 aN. N 
— = q—-KN+D|—, += — += 
.* = (i OH Gz sal 
In this equation g is taken to be the Chapman function 
z—z z—2 - 
q= aveXP| —* 7p? —cosee dexp(—=F7)| (09<¢<72) 
= 0 (7< $< 2m). (3) 
In accordance with equation (1) we shall denote by z, the height of 
maximum ion production and write, without loss of generality, 


. = 2— 2 
kK = Bexp( — H ) (4) 
where f is a constant. 
The coefficient of diffusion D may be expressed in the form b/n where 
b is a function of n and 7. For the temperatures obtaining in the F, 
region, however, the variation of b with n and T' may be ignored. We shall 
further assume that the density of neutral molecules is given by 


= Ae 
“= ng exp i } (5) 
then making use of (5), we first transform (2) by the substitutions 
Sian _2*—* 
4 exp/ i ) 


and N = vf. 


The transformed equation is 


ov . - 6 a 
en of-4..gd- sein. 
a =u "+ Gn, H? at? 
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It will be noted from equation (8) that the transformed problem is 
identical with that of heat conduction in a rod with radiation from the 
lateral surface and a heat-source distribution qf 

In equation (8) we now write 

t—kdb,  n H2/kb 
and 4(By)*l? — 2°; 


then making use of (3) and (4) we obtain 


{_ 2 cosec $) 
\ — 4(By)tf 

We require a solution of this equation which is periodic in ¢, with 
period 27. Although there are no physical boundaries, we shall see below 
that the behaviour of v at x = 0 (infinite height) must be specified. It is 
required also that NV — 0 when x > o. 


bkq, e(y'B-ix)exp (11) 


4. Solution of the equation 


The homogeneous equation corresponding to (11) is separable. Thus, 
inserting v V'(x)D(d) into the homogeneous equation we find that 
~ e—(yiBt vd (12) 
and that V satisfies the differential equation 
d?V . , 
+. (v—jfa (13) 
dx? 


where v is the real or complex separation constant. 

Solutions of the equation (13) are the parabolic functions D, s(x). Now 
the functions ; 

D(x) (n = 0, 1, 2, 3.,...) (14) 
form a complete set of orthogonal functions in the interval —« - < @. 
They are connected with the Hermite polynomials 


, a” 


H (x ye 
da" 


(e-=*) (15) 


by the equation D(x) = 2-'"e-t°H7 (x/V2). (16) 


Any arbitrary function y(x) defined in the interval —o < x2 < © and 
satisfying suitable conditions can be expanded in a series of the form 


ub(ar) S a, D(z). (17) 
n=0 


In the present case neither v nor the function on the right of (11) have 
been defined for negative values of x. However, we postulate that, at 
a great height (z: > 0, x > 0), 

O(x). (18) 
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This implies that the electron-density NV and the rate of ion-production ¢ 
are of the same order of smallness at great height. Since for all n, D,,,(x) 
is O(1) and D,,, .,(2) is O(2) when x - 0, it follows that if v can be expanded 
in a series of the form (17), only functions of odd order D,,,,(xz) can be 
present. The definition of the function zexp{—.2? cosec ¢/4(8y)!} can now 
be completed by continuation as an odd function in the interval 
2 <2 <0. Finally, the requirement that V + 0 when x - © is satis- 
fied by all D(x). 
We therefore write , 
0 = & Sansilb)Don+1(2) (19) 
Sf 


and ae a > Jon o1(%) Dy, (2) (20) 
0 


where a = (cosec ¢)/4(By)! 


] ‘ : 
= alan)! re" Dy, ,,(x) dx. 


and Jon +1(%) 
Substituting (19) and (20) into (11) we find 


> fons($)D3, (2) —(y B)* > fon ss($) Don .4(2)—} > 2 fan .3($) Don o4(2) 
n n=0 n=0 


0 


2, —$hqoely'B*) gen s1($) Dan s1(2)- 
Then, making use of the relation 
Din o4()— fx? Dy .1(2) = —(2n+4)D,, .4(2), 
we have 


x 


2 {(2n- 8) fon ss(%) ~(y B)*f on s1(4)} Don s1(2) 


n=0 
pe 2 Bh d0 (BG an +1($) Dan +12). 
Equating coefficients of D,,,,,(z) gives 
fansa(h) + (2+ 3)(B/y) fon si(b) = $hdoe(By) Gen +1()- (23) 
We now require a solution of equation (23) which is periodic in ¢ with 
period 27. The solution required is easily seen to be 


é 
fonuild) = Cc | Jon +,(x)exp{ —(2n+ 3)(B/y)*(d—a)} dat 


0 


+{exp{(4n+3)(B/y)'x}—1]-* [ gon .aladexp{—(2n+3)(B/y)(6—a)} da| 


0 


(24) 
where C = tkq,e(By)*, (25) 
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and the upper limit in the second integral is 7 on account of the definition 
(3). Hence the required periodic solution is given by 


° > feonsi($)Dan+i(2) (26) 
0 


and therefore N=C D fensil$)Donsi{2(By)*S}. 
0 


—s 
n 


5. Convergence of the solution 
It is possible by elementary methods to evaluate the integral in equation 
(22) in finite terms. Omitting the algebraic details, we find 
to n/ 9s) 3 —<nll i n 
dansa(9) = “Fi fer rs been ee) 
where Kx = (By)-*. (: 
We now examine the behaviour of f,,,.,(¢) as m > 00. We write 


( 1y"c 


2"n! 


fonsald) [Ton -1(6)+(1—exp{—(4n+ §)um})o,.1(6)], (30) 


where we have written ; (31) 


The integrals J,,,,,(¢) and J;,,,,(¢) can be expressed as Laplace aren 
¢ 
Ion 419) hy (ajer?*® da 


7 


nail?) = ( hy(x)e n{ha(x)—277} dx, 
0 


and J 


> 


2sina \3 
where h,(«) (= _ fries 


«-+sin a 


and hg(x) los(“S — *)—2u(d—a). 


K+SINn a 


It is easy to show that A,(«) can have a stationary value within the 
range of integration only if y > 8+ -' or, since is essentially positive, 
if y > 4. The values of y we contemplate using are < 4 and for such 
values h,(x) is negative and monotonic increasing. When y = 3, h(a) 
has a minimum at the end-point « = 0. The major contribution to the 
integral J,,,.,(¢), when m is large, comes from the neighbourhood of the 
upper limit. Accordingly, we expand /,(«) and h,(«) in ascending powers 
of é—« and integrate term by term. We find 


os enn _ hal) /ha(d)_, had) {hal d)}?— ha (b)ha(d)/{he(d)?? | 


n n? 


I 


2n-+ 





(36) 








DIFFUSION OF IONS IN A STATIC F, REGION 345 


where /j(¢) means the value of dh,(A—a)/da when a= 0, etc. The 
dominant term is 


F(a 


«+sind 


e- (37) 


2n 
9 gi 3 -1 

vere = eras) a) on 
The integral J,,,,,(¢) can be dealt with in a similar manner and we find 

anil) ~ €-?" "Tey 41(9)- (39) 
Now the asymptotic behaviour of D,,,,,(x), as n + 00 for fixed 2, is given 
by (6), (—1)" P(2n+2) 
ant Tin +8) 
Hence, from (24) and (40), 
fans b)Don sal) ~ 5 Fo) 4) _ tee 

x [sin(2n+ $)*x+ O(n-*)]. 

Making use of Stirling’s formula, this can be written 


°1 _ ® gin(é 3)t | 1\n+t 
fons?) Dan ia(2) ~ = Fie): aa) sin(2n-+ 3) (4) | 


(2 «+sind n n 


Dyy + 1(2) _ 


[sin(2n+ 3)ta+ O(n-*)]. (40) 





41 
where the bracketed function of n tends to unity as n > 0. vas 
It can be seen that the solution is uniformly convergent and that the 
convergence can be expected to be rapid in the neighbourhood of ¢ = 42 
(mid-day) but slow in the neighbourhood of ¢ = 0 and ¢ = a (sunrise and 
sunset). 


6. Uniqueness of the solution 

If there are two time-periodic solutions v, and v, of equation (11) which 
satisfy the conditions of vanishing at x = 0 and z = o then v = v,—», 
will satisfy the homogeneous equation corresponding to (11) and will 
likewise vanish at = 0 and x = oo. 

Now it can be seen from (12) and (13) that a purely harmonic solution 
with time-factor e'"¢ (n integral) must be of the form 

v= ree'*{AD,_,(x)+ BD,_,(—2)}, (42) 
where A and B are arbitrary complex constants, v = —i(y/8)*n, and 
D,_,(x), D,_4(—2) are independent solutions of (13). 

The condition of vanishing at x = 0 can be satisfied by taking B = —A, 
but the general asymptotic solution of D,,(—z2) (6) shows that the 
behaviour of D,_,(—2) at = 00 requires B = 0, It follows that v = 0 
is the only solution satisfying all the conditions. 

5002.47 Aa 
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USE OF GREEN’S FUNCTION IN THE SOLUTION OF 
IONOSPHERIC DIFFUSION PROBLEMS 


By J. E. C. GLIDDON 
(Queen Mary College, Mile End Road, London, E. 1) 


[Received 4 December 1958] 


SUMMARY 

The mathematical analogy between diffusion of ions in the ionosphore and heat 
conduction in a rod is used to solve the problem of diffusion with an attachment- 
type law of electron loss. The solution for the case of a constant loss coefficient is 
obtained in terms of Green’s function for the semi-infinite conducting solid. For 
the case of a loss coefficient which decreases exponentially with height a generalized 
heat-source is derived and the problem is solved in terms of the corresponding 
Green’s function. 


1. Introduction 

LN a previous paper (1), it was shown that the equation governing diffusion 
of ions in the F, region of the ionosphere was reducible to the equation of 
heat conduction in a rod with radiation taking place from the lateral 
surface. This comparison will now be used to derive solutions of the 
diffusion problem in terms of the Green’s function for the two cases when 
the loss coefficient K is (i) constant, (ii) of the form §¢*, where £ is a 
constant and ¢ is related to the space variable. 


2. Statement of the problem 
It has been shown in (1) that the electron density N(¢,¢) in the F, 
region of the ionosphere is given by 
N = lt, (1) 
where v satisfies the equation 


ev 4 ov 
a Vag 
In this equation ¢ is related to the space variable (corresponding to 
height above a fixed level), ¢ is the time measured in radians from sunrise, 
1/y is proportional to the diffusion coefficient, K is the loss coefficient of 
electrons, k is the number of seconds in a radian, and 


q = Wevte Coser} = (OS <7), 


4ykKv = —4ykql-. (2) 


= 0 (7 < ¢ < 2n), (3) 
is the rate of ion-production. 
(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3 1959) 
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It is required to find the solution of (2) which is periodic in ¢ with period 
27 and which vanishes when = 0 and when [ = o. 


3. Green’s function for the case when K is constant 
When K is constant the ‘radiation’ term —4ykKv can be removed from 


9 a 2g 2+} i » 
(2) by the substitution | (4) 


The function u(f,¢) then satisfies the heat-conduction equation 
u 


ahd C 
L(u) = —_—4y 
( (2 @d 


C7 d 


a cia 4ykql-'e kKo 


The boundary conditions to be satisfied by wu are 
u(0, dh) 0 (6) 


and lim u(l, db) 0, (7) 


Gg >~s 
while in place of the usual initial-value condition we have the condition 
of time-periodicity, 
u(l,d6+2mr7) = u(l,d) (m = 1,2, 3....). (8) 
Green’s function for this problem can be found directly by reflecting 
a unit heat source at ¢ = ¢, and time ¢ = ¢, in the boundary ¢ = 0. 
Thus, 


ss aly Weel ve). f+ E0)") 
G(l,.¢ Co» Po) a (5) Jexp{ F—e| —exp| Eee ° (9) 


The equation satisfied by G is the adjoint to (5), namely, 


2G ag 
M(G) = oat fy, ~ = 0. (10) 
Se Y 


Gveen’s theorem corresponding to (5) is then given by 


(Lu) —uMiG) déat [ [ —4ykGgl-1ekKs ddd, (11) 


where the integration is extended over the infinite rectangle 
-2n7 < ¢ < dy—e, 0<f<om. (12) 
and the limit is taken when e > 0. Integrating the left-hand side of (11) we 
obtain 
$o- € i a0 
c G — r y do-€ ba 
G pus] ddb—4y | [Gulhys de 


—2n7 7 0 
do—e ao 


(  [ —4yk@ql-ekK$ dddl. (3) 


—2n7r 0 
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The first integral on the left vanishes by virtue of (6), (7), and (5) 
letting « > 0 we obtain 


u(Los $o)— [ @(L, —2nm | Lo, bo)t(Lo, —2nz) aL 
0 


do a 
=k | &k$ dg J GUL, $ | Los bo)ah-t a+ 


2a 


+ 5 kf Akis-2m ag i G(L,6—2mm| Ly. dq)ql-* dl. (14) 
m=1 6 


Here we have made use of the periodicity of g and it may be noted that 
the limits (0,27) with respect to ¢ can be replaced by (0,7) according 
to (3). 

We now evaluate the integrals on the right-hand side of (14). We have 


k [ G(Z,4\% -1 dt = kq,en-* Y } (5. a) 
Haein i: Nae P\So—4) * 


0 


x [Sex| —<*(5% + coseed) + F ye 5st} 


2y¢ 
tex {- eases )—7r, i] dt. (15 
Wad” a ceses tare ” 
Integration can be effected with the help of the error-function integrals 
giving 


[tea dt— | tea dt = 2 seme lertol —5-)+ erfe(-)} 
a d 2a \ 2a 


b 
= gt 


eb?/4a? 
2a* 


Applying t!.s result to the right-hand side of (15) we have 
3 haa 3 
k [ GL.4 \Lododae-* de = bgyet,(1-+2°—* cose) * x 
0 
(2 cosec b 
1+[(¢o—¢)/y]eosec ¢ 
Similar results are obtained for the remaining integrals on the right of 


(14) by writing ¢—2mz in place of ¢. Also, letting n > o and using the 
periodicity of u(¢,¢), we see that the integral on the left of (14) tends 


x exp} — 





(17) 
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uniformly to zero. Finally, interchanging 5, ¢) and ¢, ¢ we obtain the 
required result, 


é , 
u(l,d) = kqge€ | Aj "exp[EK4, — > cosec $a} dd + 
a ( 


C?cosec do) 
aaa | 


n 


A, exp| kK ($y—2n7) 


ny 
0 


dd, (18) 


=1 


> 

. q aN 

where A 1, ro —— cosecd, (n = 0,1, 2,...). (19) 
7 
/ 


The electron density V, as given by (1), (4), and (18), agrees with the 
result previously obtained by Ferraro and Ozdogan (2). 


4. Green’s function for the case when kK = B@? 
When the loss coefficient decreases exponentially with the height we 


write (1), kK = pt, (20) 


where § is constant. Equation (2) is then transformed by the substitutions 


4(By)tl? = 2°, 
and 


Writing, for convenience, 


we obtain 


The differential equation to be satisfied by the Green’s function G is 


8(t)—7)8(x—a9). (26) 


eG 
oe 
ox” 


In the first place we require the solution of this equation, in the range 
0 <x < ®, which vanishes at x = --oo. Assuming an expansion for 
G in terms of parabolic functions D,(x), we write 


co 


G(x, T | Xq; To) 2, Pnlt) Dal ). 


Substituting this series into (26) we find, after simplification, 


> (Pals) —(N+4)pp(7)}-Dy(x) = 8(t>—7)3(x—29). 


n= 
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Multiplying both sides of (28) by D,(x) and integrating from —o to 00 we 


find that 
, ; D, 
P(t) —(N+ 4) Pa(7) = 8(t9—7) mer 


e-n+tXre-7) D (a9) (30) 


(29) 


whence p,(t) = (nin! 
and so, formally, 
is) 


. ] 1 a, 
Ga, 7 | qT) = ay aa iXte-7) D) (275) D,, (x). (31) 
Again G(x,r | 2,79) can be expressed in closed form by use of the series 


of products of Hermite polynomials, viz. 


(32) 





> G2 H(e)Hy) = (1-24) exp} 
L— | n. 
From (31) and (32), with the aid of the relation 


D, (x) = 2-*"e-** H, (x/~2), (33) 


Qayz— (x?+-y?)a* 
1—z2? 


we find that 
(x?+-22)cosh(7+,—7) — 2xay 
4sinh(r,—7r) 





G(x,7|%,T)) = {4asinh(t,—7)}-* exp| — 


| (34) 


It is easily verified that 
@ 


lim | G(x, 7 | 2,7) dz = 1, (35) 


ites 

so that (34) is the analogue of the unit heat source used in Section 3. 
Also, as in the previous section, the Green’s function which vanishes 
when x = 0 is given by reflecting (34) in x = 0, that is, returning now 
to the use of ¢, do, 

G(x, p | Xo, bo) — G(x, | —Xq; do). (36) 
The solution of (25) is now obtained as in section 3, Green’s theorem being 
of precisely similar form. We find that 


do @ 
(29 bo) = HhqoeBtytlim [pw dg | xe-testooneed 
—~@ 0 


—2n7 


x {G(x, db | Xo, $9)— G(x, b | —%o, bo)} dx, 


or 
U(Xo, do) 
- tkqge (By)-* | d¢ | ae~ixx* cosee $f (ar, ¢ | Xo, oo) — G(x, ¢ | —Xo, $o)} dz+- 


0 0 


+ Hho e(By) > 


x S| [ xe tertooneed (G(x, —2ne | x9, $o)— G(x, d—2me | a, $o)} de. 
00 (37) 
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We find with the help of (16) that 


x 


re erent Gis, ¢ Xo; Py) — G(x, —Xo,$o)} dx 


- 4 Se x5 T'(do, $) 
XA (9, d)exp P\— 4 (do, $) 


(bo, ¢) = « cosec ¢ cosh u(d,—¢)+sinh p(}d,— 4), (39) 
A(¢o, ¢) = « cosec ¢ sinh p(¢d,—¢)+ cosh p(d,—¢). (40) 
In a similar way we can express the integrals under the summation 
sign in (37) in the form 
P { «2 T'(bo,¢—2nz) 
| Xi A(do, d 2nz7)} texp)—~° 4 A a ae rt d 
0 


Finally, interchanging x,, ¢, and 2, 6 we have 


j | x T'(¢, bo) | 
v(x, b) = tkqge(By)-* { x{A($,4,)}texp| —— ES! dy t 
2 10 py) ! 0 P 4 (, bo)I 0 


(41) 
d 


7 


2 I'(¢, dy — 2nz) Na 
Llhkq,¢ y> f —2nz)\-! = 42 
bkqy e(By) 2 J a{A(d,bdo—2n exp) —F N6.dq—2n7)| ddy. (42) 


5. Proof of equivalence of the two solutions obtained for the case 
kK = pt 
We now show that the solution obtained in section 4 is equivalent to 
the solution obtained previously in (1). 
The series form (31) of the Green’s function shows that 
G(x, | %,b9)— G(x, d | —Xpq, bo) 
a, 
ee bee ak n+ Bylge— —O) Dan (2% ») Dania (2). (43) 
(277)? int, (2n+1)! 
n= 


Using this form in (37) we have 


: B.. " 
U(Xp, ho) $kqye(By) i S < . = are! | « ean +Duide-) dd > 
(27) (2 n+l)! | 
0 


ben 


‘| qge-tKr? coneed De. . .(z) dx- khqy¢ (By)-* S 2 i 


n=08 


yt (2 n+1)! 
" ° 

x | e-2n+imdo-d+297) dd ze txx* cosec p Dry 44(") dx. (44) 
0 0 
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In the notation of (1) we have 


| ae-ixz*coseed J) (x) da, 


0 


2 l 


Jonsi(o) = (Qn) (Qn+1)! 


so that (44) may be written 
to 
U(X, bo) = $kqye(By)* 2, Dyn, 44(Xpq) ( Janse +Dude-#) db +- 


0 
«o 3] by 
+ tkgqge(By)* > > Dan sao) [ Gan+a(p)e-2n+Pmide-d+20m dd. (46) 
n=08=1 0 
tearranging and interchanging 2, ¢, with x», ¢), we have 


o{? 
v(x, ) - kkqye(By)-* 3 | | Jon+1( Poe 2" Hoo) dd.+ 
n= 6 


v 
7 > {e“an+Bym}s go 1s (pyle an to-do) Abo} Dyn (2) 
é s=1 


a 
or, inserting the value of > {e-@"+9#7}s, 
s=1 


o(? 
v(z,¢) = $kq,e(By)-* 2 {Foam sal de-nstmd-$0 ddy+ 
n=0\9 


+ [ee +8e" — 3) -2 | Jan+1(Pole~2" + Deo —-oe) a Dyn i3(%). (47) 
0 
This agrees with the result obtained in (1) when written in the form 


v(a. 4) = ¥ fansi(@)Dan al). (48) 
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SUMMARY 


Some inequalities relating to magnetism are proved and their relevance to the 
interpretation of local anomalies in the earth’s magnetic field is discussed. From 
observations of the resolved part, in any one direction, of the earth’s magnetic field, 
these inequalities could be used to calculate (subject to certain hypotheses) a maxi- 
mum possible value for the depth below the earth’s surface of the top surface of 
the magnetized body producing the anomaly. 


1. Introduction 

In this paper some inequalities are derived which may be of use in 
discussing the results of surveys of local anomalies in the earth's magnetic 
field. The problem of interpretation, which is to determine the nature 
and structure of the magnetized bodies causing the anomaly, is complicated 
by the fact that the distribution of dipoles which will produce a given 


field at the earth’s surface is not unique. The methods normally used in 
discussing this problem (see (1), (2), (3)) draw on information from other 
sources, often geological, to form a rough model and then adjust the 
characteristics of this model to fit the observed anomaly. When there is 


only little information available from other sources not much can be said 
with confidence about the body causing the anomaly. It is in such situa- 
tions as this that the method described below might be of use. It produces 
some useful information about the body but requires only light assump- 
tions to be made about it. The method can be applied to the readings of 
surveys using either a vertical magnetometer or a total field (airborne) 
magnetometer. The analogous problem for gravity anomalies has been 
discussed by Bullard and Cooper (4) and by Bott and Smith (5) but neither 
of these methods can be used for the study of magnetic anomalies. 

The writer is indebted to Dr. M. H. P. Bott for suggesting the problem 
to him and for much helpful advice on various aspects of geophysical 
prospecting. 

2. Summary and discussion of the results 

Take a system of rectangular cartesian axes Oxyz with the y-axis 

pointing vertically upwards and the other axes in the plane of the earth’s 


{Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959] 
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surface. Let v be any fixed vector and let A(x, z) denote the resolved part 
parallel to v of the magnetic field produced at the point (z, 0,z) by a finite 
magnetized body B which lies beneath the earth’s surface. The body B 
will be said to have parallel magnetization if the intensity of magnetization 
I has the same directiont and sense at every point of B though its magni- 
tude |I| may vary from point to point. We prove the following results: 


THEOREM |. Jf the finite body B has parallel magnetization and lies wholly 
underneath the plane y = —h(h > 0) then 


hf [f(e)Pdx <4 | [FP de, 


—@ —@ 
2 


+ 
where f(x) = f A(x, z) dz, 


—@® 


and F(x) = [ f) dt. 


Also AV(F2) < 28 | [ F(x)}? de, 
where V(F?) denotes the total variation of the function | F(x) |* over the range 
-O < 2% < +0. 


The number V(F?) is easily calculated from a graph or table of the 
function by means of the relation V(F?) = 2(M—m), where M denotes 
the sum of all the maximum values of | F(x)|* and m denotes the sum of 
all its minimum values. 


THEOREM 2. If the finite body B has parallel magnetization and lies under- 
neath all three of the planes y = —h, y = 34x, y = —3+z, then 
h? ro dx< 3 [ f(x)? dx, (2.5) 


+o 

and hV(f2) < 34.2 f (f(x)? dx, (2.6) 
where f(x) is given by (2.2) and V(f*) denotes the total variation of the 
function | f(x) |* over the range —oo <4 < +00. 

The number V(f?) can be evaluated in the manner described above for 
V (F?). 

If the vector v is taken to point vertically downwards, then A(z, z) is 
the magnetic anomaly at the point (x, 0,z) due to B, as obtained from the 


+ The vector I is not required to have the same direction as the vector v of course. 
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readings of a vertical magnetometer. If the magnetic ground field of the 
earth can be regarded as constant in both magnitude and direction over 
the region considered and the vector v is taken parallel to this ground 
field, then, to the first order, A(x, z) is the magnetic anomaly at the point 
(x, 0,z) as obtained from the readings of a total field magnetometer pro- 
vided that this anomaly is small. If the function A(zx,z) is completely 
known, the integrals in (2.1), (2.4), (2.5), and (2.6) can all be evaluated. 
Any one of these inequalities then provides a maximum possible value 
for the depth below the earth’s surface of the top suriace of the body B. 
The actual depth A can be any value less than or equal to this maximum. 
It should be noticed that in proving the theorems no assumptions are 
made about the shape of the body. The inequalities (2.1) and (2.5) are 
‘best possible’ in the sense that the equality sign is actually attained when 
the body B degenerates into a single dipole lying in the plane y = —h. 
Our method of proof suggests that (2.4) and (2.6) are not best possible 
in this sense. For the case when the body B is a single dipole it can be 
shown that they overestimate its depth by 12 per cent and 33 per cent 
respectively. The hypotheses of Theorem 2 require that the body B lie 
underneath a wedge of angle 60° formed by the planes y = +-3!x. This 
is a serious restriction on the usefulness of the result. It is not clear from 
our method of proof that this is the bluntest wedge under which the 
inequalities (2.5) and (2.6) are valid but it can be shown that (2.5) is not 
valid if the wedge angle is increased to 82°. 

A body B* will be described as two-dimensional and extended parallel 
to the z-axis if it can be regarded as a bundle of uniform, linear, dipoles 
each of which lies parallel to the z-axis and the body extends to infinity 
in both directions. Both the intensity of magnetization I inside such a 
body and the magnetic field H produced by it, are functions of x and y 
only and are everywhere perpendicular to the z-axis. The magnetic 
anomaly A(x) which is the resolved part parallel to v of the field H at 
the point (x, 0,z), is clearly a function of x only. Theorems | and 2 cannot 
be applied as they stand to two-dimensional bodies since the integral in 
(2.2) is then always divergent. For such bodies the following result will 
be shown to hold: 


THEOREM 3. If the body B* has parallel magnetization, lies underneath 
the plane y = —h, and is a two-dimensional body extended parallel to the 
z-axis, then the inequalities (2.1) and (2.4) are both true if f(a) is taken to 
mean A(x),and F(x) is given by (2.3). If, in addition, B* lies underneath 
both the planes y = +-3'x, then the inequalities (2.5) and (2.6) are also true 
with this meaning for f(x). 
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All these results assume that the body producing the anomaly has 
parallel magnetization. Whilst this will never be satisfied exactly, experi- 
mental evidence (see (6)) indicates that it is a fairly good approximation 
to the state of magnetization in most igneous rock masses. In proving 
Theorem I, no assumption is made about the orientation of the coordinate 
axes in the plane of the earth’s surface. In practice, however, it is 
desirable that the z-axis should not be taken parallel or nearly parallel to 
the direction of magnetization I. As explained at the end of section 5, 
any deviations from the hypothesis of parallel magnetization can be ex- 
pected to produce a more marked effect on the results when I is nearly 
parallel to the z-axis. In general I will be inclined to the horizontal plane 
at an appreciable angle and then no caution will be necessary in choosing 
the direction of the axes in the earth’s surface. 

The procedure in this paper will be to prove Theorem 3 first and then 
to deduce Theorems | and 2 by a simple extension of it. The key to this 
extension is the following result: 


THEOREM 4. Let H(x,y,z) be the magnetic field produced at the point 
(x,y,z) by a simple dipole situated at the point (£,n,¢), whose moment is the 
vector M = (A,y,v). If we write 


H*(x,y) = | H(x, y, z) dz. (2.7) 


then H* is equal to the magnetic field produced at the point (x,y,z) by a uni- 
form linear dipole which passes through the point (€,,¢), lies parallel to the 
z-axis, extends to infinity in both directions, and has for its moment per wnit 
length the vector M* = (A, p, 0). 


It should be noticed how this result is used in section 5 to extend 
Theorem 3. Any other method which may be found for estimating the 
depth of two-dimensional bodies from observations of the magentic ano- 
malies produced by them, can clearly be extended in the same way to 
produce a method applicable to three-dimensional bodies. 


3. Preliminary calculations 


In this section some preliminary results are proved. 
Lema |. Jf X and ¥ are real and X > 2h > 0, then 


(2h)-*re(X+iY)— > re(X+iY)-. 
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Proof. If tan@ = Y/X, then X+iY = Xe'®/cos6. Hence, 

(2h)-*re(X +iY)-1—re(X +-iY)-3 = (2h)-? re(X-1 cos 6 e-!?) — 

—re(X-3 cos*# e-3*9) 
ay ! — (F) cos @ cos 30] 
> 0. 
Lemma 2. If X and Y are real, X > 2h > O and |Y| < Xtan}nz, then 
(2h)-? re(X+i¥)-3 > re(X+1tY)~. 

Proof. If Y/X = tan@, it may be supposed that |@| < 47 since 
Y|/X <tanj}r. When this is so, cos3@ > 0 and cos3@ > cos5@. The 
proof of Lemma 2 now proceeds like that of Lemma 1 and it is omitted 
here. 


Lemma 3. Let d(z) = S(z)/T(z), where S(z) and T(z) are polynomials in 
z and the degree of T(z) exceeds that of S(z) by 2 or more. If $(z) is analytic 
at all points of the real axis, then 


+ @ 


( d(x) dx = 2zip, 
where the integral is taken along the real axis and p denotes the sum of the 
residues of $(z) at its poles in the half-plane imz > 0. If there are no poles 
in this half-plane, then p = 0. 


This result is a special case of a well-known theorem in the calculus of 
residues. 


Lemma 4. Let z = x+iy and let a,, dg,..., Ay, 5y, bg,..., 6, be complex 
constants with ima, < 0 forallr. If 
=. Ob 
Qz)= 5 —, (3.1) 
S- z—a, 
TP |dQl? a. bb 
then |\——~| dz = 2n(2p)! > -———_—__*___.., 3.2 
. . | dz? wv? = 2 [ta,—14, |*? +? ; 


where the integral is taken along the real axis and G,, b, denote the complex 
conjugates of a,, b, respectively. 


Equation (3.2) is true for p = 0 provided that in this case d?Q/dz? is 
taken to mean Q and 0! = 1. 
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Proof. If the bar denotes the complex conjugate then 


“Qo d°Q\ (0"Q 
i wi - | (GS\G3\ 


y (—1)?p!b,\{ Ss (—1)>p!b, 
ap 2 eo a,)?*| |, @4r ee 


+o 


SF (5, [ [@—a)e—a)} dx. (3.3) 
r=1 s=1 


—-@ 


I 


Now |(z—a,)(z—d,)|-?-! is a rational function of z which is analytic on 
the real axis and everywhere in the upper half plane except at the point 
z = a, where it has a pole of order p+1. Using Lemma 3 it can be shown 
that re 

[ [(x—a,)(e—4,)]} #1 dx = 2n(2p)!(p!)-*(ia,—ia,) 2? 
Equation (3.2) is now obtained by substituting this value of the integral 
into (3.3). j 


Lemma 5. Let the function Q(z) be as described in Lemma 4 and let 
d?Q/dz” = U+iV, where U and V are real. Then 


[ va — =} ji ora 


—-@o _- 





where the integrals are taken along the real axis. 
Proof. If $(z) = {d°Q/dz?}* = (U+iV)* 


then ¢(z) satisfies all the conditions of Lemma 3 and has no poles in the 
half-plane imz > 0. It follows from Lemma 3 that 


ak 
= | de) de 
+0 +0 

m f (U2—V?2) da+i f 2UV dz. 


It now follows from the real part of this equation that 


{ Utd “= [7 dx = +f (UP dx. 
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4. Two-dimensional bodies 
In this section Theorem 3 will be proved with the help of the following 
result: 


Lemma 6. Consider a collection of n uniform, infinite, linear dipoles with 
moments M,, M,,..., M,, respectively, which are all extended parallel to the 


rT 
¥ 


f-axis of a system of rectangular axes Ofnl and cut the plane ( = 0 in 
the points (£,,7,,9), (&, m2, 9),..., (€., 2» 0) respectively. Suppose further 
that », is negative for all r and that all the moment vectors of these dipoles 
point in the same direction, which is perpendicular to the C-axis and makes 
an angle « with the direction of the positive €-axis. Let v be a fixed vector 
whose direction cosines with respect to these axes are (A,,v), and let A(x) 
denote the resolved part parallel to v of the magnetic field produced at the 


zr 
point (x,0,0) by these linear dipoles. If F(x) = { A(x) dx, then 
l — (d? F\? n n . 
nr Mod = > ¥MM. —ié,—7,—7,|~2?-1. 
47(2p)!(A?+-p?) \ dx? | a r x | r .re[ig, if, Vr ns] 
x (4.1) 


. 


This is valid for p = 0 provided that in this case d? F/dx” is taken to mean 
F and 0! ==s 2. 


Proof. The two-dimensional magnetic field produced by this collection 
of linear dipoles has a complex potential (see (7), p. 182) whose value at 
the point (zx, y, 0) is 


=. 2M, e% 
W(z) = p3 .<, (4.2) 
r=1 we 
where z = x+y and a, = é,+in,. In terms of this complex potential, 
the resolved parts parallel to the coordinate axes of the field at the point 
\*, 9, 0) ane (—redW /dz, imdW/dz, 0). 


The resolved part parallel to v of the field at the point (x, 0, 0) is therefore 


A(x) Are tid +pim iad = re{C ry 


where C = —(A+iu). It follows that 


>| dx = re{CW(2)}, 


since W(x) >0 as x > —oo. If b, = 2CM,e, then CW(z) as given by 
(4.2) is identical with the function Q(z) defined by (3.1) and so 


F(x) = [ re(C 


aw) 
U 


F(x) =reQ(x) and d?F/dx? = re(d’Q/dz’). 
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From Lemmas 4 and 5 it now follows that 

7 2 ~~. 2 

ig is} ¢@ v=3( ole 
2 dz? 

x o“* 


~ore > > aa —1a,| 1 


r=1 8s=1 








On substituting for a, and 6, in this and remembering that YM, is real, we 
obtain 





(d?F\? , aciQie a 
\ dx | Caer > > [1€,—1€,— —1, jr) 


r=1 s=1 


The result (4.1) follows on taking the real part of this equation. 


Proof of Theorem 3. We shall prove Theorem 3 first for the case when 
the two-dimensional body B* consists of the collection of n linear dipoles 
described in Lemma 6. If these linear dipoles all lie underneath the plane 
n —h, then Ne S —hforallr. If 


1£,—16,—N— Ns a X+iY, (4.3) 
then X = —n,—7», > 2h > 0. Lemma | now shows that 
(2h)-* re(ié,—i€,—,—,)71 > re(ié,—i€,—7,—7,), 


for all r and s. It follows that 
(2h) > > M, M, re(it,—1£,—4,—4)* 
= >> > M, M, re(igé,—t£,—,— ”)*. 


If we substitute for the two sums from (4.1) this becomes 


+@ 1 ari 
ony-2 ( F 1 (4F\* a ; 
(2h) [ Pae>s | iz| @ (4.4) 


Since f(x) = dF /dz, this is the inequality (2.1). 
The inequality (2.4) is now deduced from this. The total variation of 
F? is given by 


+0 
V(F%) = | ae de 








@ 
5092.47 Bb 
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Schwarz’s inequality (see (8), p. 381), when applied to this integral, shows 
that 


+@ 


ee ee 
[V(F?)] ‘| I dx \dz| 


a> @ —@ 


dx. 


Substituting into the right-hand side from (4.4) we obtain 


+x 


[V(F?)/? < 2h | ( F? ae 


a 
The square root of this is (2.4). 
It will now be assumed that the n linear dipoles described in Lemma 6 


all lie underneath both the planes 7» = +3#€. Then, in addition to the 


inequality », << —h, we have —», > 3*\€,| for all r. If X+7Y is given 


by (4.3) then as before X > 2h > 0 and also 
xX Np—Ne > BHE,|+-34/é,| > 3tE,—E,| = 3t/¥ 
which is equivalent to Y| < Xtan47. Lemma 2 now shows that 
(2h)-* re(ié, —i£,— n,— ,)-3 > re(ié,—t€,— ,— 9,)~° 
for all r and s. It follows that 
(2h)-2 SSM, M, re(ié,—i,—9,—1,) 
> > Dd M, M,re(té,—t£,—,— 1). 


On substitution for the two sums from (4.1) this becomes 


(2h) 


“(APY a, | i | Gal 


2! | \da | \ dx? | 

Since f(x) = dF /dx, this is the inequality (2.5). Finally, the inequality 
2.6) is deduced from (2.5) in precisely the same way that (2.4) was 
obtained from (2.1). 

Theorem 3 has now been proved for the case when the two-dimensional 
body B* consists of a finite collection of linear dipoles as described in 
Lemma 6. A continuous two-dimensional body B* can be subdivided into 
a large number of infinite magnetized rods which approximate to linear 
dipoles. The magnetic field of B* can therefore be imitated as closely as 
we please by the field of a finite collection of linear dipoles for which the 
inequalities of Theorem 3 hold. These inequalities are therefore also true 
for the field of the continuous body B*. This intuitive argument can be 
justified easily by appeal to the general convergence theorem of Lebesgue 
(see (8), p. 345). This completes the proof of Theorem 3. 
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5. Three-dimensional bodies 
In this section Theorem 4 is proved and then used to extend Theorem 3. 
Proof of Theorem 4. Let i, j, k denote unit vectors parallel to the 
coordinate axes. If r = (x,y,z) and a = (€,,¢) then the magnetic field 
H of the simple dipole M is given by H = —VQ where 
Q = —M.V\r—a|- 
(see (7), p. 133). From (2.7), 


H*(z,y) = — i. ae. +4 S)oe, y,z) dz 


ox cy 


a -(i a+ix) | Q(x, y, z) dz 
ox * Oy 
= —VO*, 


~% Q(a, y, z) dz 


where 
| (x —~ + p- “+¥2)ir— a|-! dz 
; cy oz 


Since the integrand in the last integral contains the variables z and { as 
a function of z—{ only, the integral can be rewritten as follows: 


Q*(x,y) = — | ( z +5 )ir—a)- dt 
Cx cy 


[ M*.V\r—a|-1 dz, (5.2) 
—@ 
where M* = (A,y,0). Equation (5.2) shows that * is the potential of 
a linear dipole which lies parallel to the z-axis and has moment M* per 
unit length. Equation (5.1) then shows that H* is the field of this linear 
dipole. This completes the proof of Theorem 4. 

Proof of Theorems 1 and 2. Let H(x,0,z) be the total field produced 
at the point (x, 0,z) by a collection of n simple dipoles with moments M,, 
M.,,..., M,, respectively, which all lie underneath the plane y = —h. If 
H*(x) is obtained from H(z,0,z) as in (2.7) then Theorem 4 shows that 
H*(zx) is the field produced at the point (x, 0,z) by a family of linear dipoles 
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whose moments per unit length are Mf, M$.,..., M* respectively. These 
linear dipoles all lie parallel to the z-axis and underneath the plane 
y = —h and if the vectors M,,..., M,, are all parallel, the same will be true 
of the vectors Mf...., M*. When this is so they form a two-dimensional 
body of the type envisaged in Theorem 3 and if f(z) denotes the resolved 
part of H*(z) parallel to the fixed vector v, then f(x) satisfies the in- 
equalities (2.1) and (2.4). If it is assumed in addition that the nm simple 
dipoles M, M,, all lie underneath the planes y = + 3!z, then the linear 
dipoles Mf,..., M% also satisfy this hypothesis and the function f(x) satis- 
fies the additional inequalities (2.5) and (2.6). The resolved part of 
equation (2.7) parallel to v is f(x) = { A(x,z)dz, where A(x, z) denotes the 
resolved part of H(z, 0,z) parallel to v. This shows that f(x) is the function 
defined by (2.2). Theorems 1 and 2 have now been proved for the case 
when the body B consists of a finite collection of simple dipoles all pointing 
in the same direction. The results can be extended to the case when B 
is a finite continuous body with parallel magnetization by the same 
argument as that used in the proof of Theorem 3. 


6. Nearly parallel magnetization 
Suppose now that the moments M, M,, of the » simple dipoles 


referred to in the previous proof are not exactly parallel to each other but 
are very nearly so. Theorem 4 shows that each vector M* is the ortho- 
gonal projection on the plane z = 0 of the corresponding vector M,. If 
the angle between M, and M, is small then the angle between their ortho- 
gonal projections M*¥ and Mf will also be small, except in the case when 
M, and M, are nearly parallel to the z-axis. If this possibility is excluded 
then the moments Mf,..., M* of the corresponding linear dipoles will also 
be very nearly parallel. For this reason it is important that in practical 
applications the z-axis should not be taken nearly parallel to the direction 
of magnetization. 
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SUMMARY 

In this paper two iterative methods are developed for the numerical solution of 
polynomial equations with real coefficients. The first (the 2) method is based upon 
a systematic application of Routh’s test (1) for the number of zeros of a polynomial 
with positive real parts. The second (the a) method uses the test in a different manner 
to evaluate quadratic factors directly. The methods have certain desirable features: 
(i) they always converge, even from an arbitrary start, (ii) numerical checks can be 
used at each stage, and (iii) with modern fully automatic desk machines with retain- 
ing keys they are relatively quick. Some properties of Routh’s test are proved and 
an algebraic proof of the test itself (and two proofs of the fundamental theorem of 
algebra) are also given. 


1. Introduction 

Tue methods given in this paper are partly developed from the test func- 
tion methods of Frazer and Duncan (2). They have the advantage over 
the latter of not requiring any preliminary investigation by trial and error 
to obtain a near approximation before converging rapidly; in fact they 
always converge, even from an arbitrary start, and so have also this 
advantage over methods such as Lin’s (3) and Friedman’s (4) which may 
fail to converge however close the first approximation used to start. Lin’s 
method has been modified so as to give greater probability of convergence 
by Aitken (5). In the later stages of the application of the a and x methods 
inverse interpolation from a suitable test function enables rapid con- 
vergence to be achieved, especially if Aitken’s (6) or Neville’s (7) repeated 
linear interpolation form of Lagrangian interpolation is used. 

In section 2 the x method is first described. This applies Routh’s test 
directly, with suitable changes of origin, and determines first whether real 
and/or complex zeros are present. If real roots are shown to exist they 
are then found in the usual way, while complex roots are found from the 
calculations involved in applying Routh’s test itself. 

The a method is described in section 3 and owes its origin partly to 
a method for solving quartics given by Mack and Porter (8). It consists in 
finding quadratic factors of the polynomial directly and applying Routh’s 
test in an indirect, but equally powerful, way to do this systematically. 

(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959) 
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Complex zeros are found by the a method with appreciably less computa- 


tion than by the x method, while real zeros can be found by one extra 
calculation at each stage (or can be found in pairs as quadratic factors). 
For degrees of six and upward it is usually the more economical method 
particularly if all zeros are known to have real parts of the same sign, or 
if the first application of Routh’s test shows this to be so. 

The two methods may be combined in a number of ways. Thus one or 
two applications of the x method may be used to convert the equation 
into one whose real parts are all of one sign, whereupon the a method can 
be used to complete the solution. It will also be shown that some of the 
quantities needed in the a method are computed incidentally during the 
development of the x method, and time and labour may be saved by using 
them. The methods also combine well with Bairstow’s (9) method of 
improving a quadratic factor which usually needs a three or four significant 
figure approximation to start with but then converges rapidly. 

At present the method recommended by many computers is the root- 
squaring process due to Graeffe (10). This is difficult to use without 
making mistakes; as significant figures are usually lost, and complications 
arise if there are several complex zeros present though Olver (11) does 
give a technique which avoids the main cause of error. Such difficulties 
are avoided with the methods of the present paper. 

The time required by a particular method is an important criterion of 
its usefulness and some average times are given in the text. The methods 
can also be programmed easily for calculation on electronic computers. 


2. The x method 
2.1. Routh’s test 
We begin by stating Routh’s test (1). Suppose that 
f(z) = coz +c, 2*-!+-¢,2"-*+...1¢,,- ; (1) 


is a polynomial in z of degree n with real coefficients. Routh’s test consists , 
essentially, in taking the two polynomials formed by the even and odd 
powers of f(z), dividing that of higher degree (i.e. cyz"-+-c,2"-*+... = fy, 
say) by the other (c, z”-!+-c,2"-3+-... = f,, say), obtaining a remainder f,, 
then dividing f, by f., which is of degree n—2, obtaining a remainder fs, 
and so on until f,,, a constant, in fact equal toc,, as we shall see, is obtained. 
This h.c.f. process is analogous to the Sturm process for finding the number 
of real zeros. Then, as Routh (1) showed, the number of changes of sign 
in the coefficients of the leading (i.e. the highest power) term in each f; 
equals the number of zeros of f(z) with positive real parts (a proof is given 
in subsection 4). 
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In the form of detached coefficients this process yields the array 
C4 C6 . . » 34 
C5 
cy c®) 


cf?) 


(n—3) (n—3) 
Cn-3 Cn—-1 
(n—2) (n—2) 
Cn-2 Ch 
(n—1) 
Cn-1 


(n) 
Cn 


— . . 0 am gli) (r) i i 
where, with c., = ci’, co,,, = ci)... the ci”),, satisfy the relation 
7 ae —2) (r—2) /o(r—Wielr—1) 

Cr+2s aps Cr as —\Cr—3 Cry Jeeta (3) 
provided r+-2s <n—1. In the final column of (2), where r+-2s = n, it 
can be seen that all terms are the same, that is 

A) =z COW ae COW® a .., = C,,. (4) 
[t is worth noting that each row of (2) is generated from the two rows 
immediately above it in the same way and, also, that labour can be saved 
in computing the rth row by first calculating c{=?/c"> and then retain- 
ing it as a multiplier. 
The number of changes of sign in the sequence 
Coy Cop CO, CP ,.005 CO, 200y CP (5) 
then gives the number of zeros of f(z) with positive real parts. 


A check on the accuracy of the calculation of the array (2) may be made 
from the following result: If s” is the sum of the rth row, i.e. if 


a = cl) +c) +c) 44+... (6) 
then it can be seen from (3) that 


= —-2 (r—2) /p(r—1 -1 
GH) == glt-2 fel? = 2) [ol “Mh gr -D, (7) 


2.2. Description of the x method 


In the x method we change the origin from z = 0 to z = x, compute 
the new coefficients of f(z) when expressed in powers of (z—zx), and then 
apply Routh’s test to these coefficients. Thus, we first find the coefficients 
C;, where 

o 


f(z) = Oy)(z—2)" +C,(z—x)"-1+C,(z—2)*-? +... 4+-C,_,(z—2)+C,. (8) 


This can be done by repeated division of f(z) by z—z. An alternative, and 
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quicker, method when retaining keys are available, is to compute them 
by calculating column by column the following array: 


(¢, "10 (Co =)og (Cg =May - + + (Cn-y=)ln-r9 (Cy =)eno 


Coy 3, be 1,1 Cc C,,) 


Co 


2 
e 


“n—1,21 


(10) 


In words, each term in (9) is the sum of the one above and the product 
of the term to the left times x. Accuracy may be checked by using the 
relation ’ 

+0, = f(x+1) (11) 
which follows from (8) when z = «+1. In practice, it is convenient to 


sum ©, |-C,+C,+... and C,+C,+... separately as these sums are useful 


. 
in checking the Routh’s test calculation which follows. 

Routh’s test (see (2)) is now applied to the C;. Denote the general term 
in the consequent array (2) by C\’),,. Then the number of changes of sign, 
pix} say, in the sequence C!” (r = 0,1, 2,...,”), gives the number of zeros 
of f(z) with real parts greater than x. 


Once we have found 2,, x, such that p{x,} + p{x,} then we know that 
at least one real part, x) say, lies in the range [z,,x,], and further sub- 
division of this range will find x, as accurately as desired. We can, in fact, 
specify such a range without trial, for if z, is the zero of f(z) with least 
modulus then Z| < |e,/eg)" = X, say, (12) 
and a real part must therefore lie in {—X, X]. Sometimes too it is known, 
or the first test (which is usually with z = 0) shows, that p{0} = 0 or n. 
In both these cases all real parts must have the same sign and, hence, one 
real part must lie between 0 and —c,/(nc,), a range which may be narrower 
than [0, X] or [—X, 0]. 

A given range [x,, x,| can be subdivided in half and Routh’s test applied 
again, and so on. If no indication that the true real part x, is nearer one 
boundary than the other, this subdivision by half can be repeated a few 
times. Then it is usually possible by inverse interpolation from a suitable 
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test function to find x, with great rapidity of convergence. This we now 
discuss. 

The procedure differs in the two cases of real and complex zeros. We 
consider these two cases separately. If p{z,} and p{x,} differ by an odd 
number it will be found (see subsection 4.1) that f(z,) and f(x,) have 
opposite signs. This is a well-known condition for a real zero to exist in 
[2,,2,| and it can be found in the usual way. In the later stages of this 
process inverse interpolation from the values of f(x) itself can be used. 

If p{x,} and p{x,} differ by an even number some suitable point x, in the 
range [x,,x2,| is chosen and p{z,} found. Subdivision is carried on until 
a range is found in which either p{x} changes by an odd number or p{z} 
changes by two, whilst f(z) remains unchanged in sign (or appears to do 
so). If f(x) does not change sign in the interval then it will be found that 
at one and only one intermediate point x), C\">" (the penultimate quantity 
in the Routh test array (2)) will vanish. At this point it will also be 
found that C1n-2/C(M=2) — CM=9/CM=9 = yZ,_ say, (13) 
where y2 is positive. Finally, it can be said that x,+iy, are zeros of f(z). 
It should be noted that it is possible for f(z) to have two real zeros close 
together and that it is then sometimes difficult to be sure of the nature 
of the zeros until somewhat later in the proceedings. In such cases 
C™ = f(x) is usually rather small and this gives an indication of the 
possible existence of two close real zeros. Otherwise it is generally safe 
in the later stages to interpolate inversely from C(">" to determine 2p. 

Sometimes there is a loss of significant figures in C(*>) and then accuracy 
is gained by calculating the last three rows of the test array (2) as follows. 
Evaluate CM=9/CM=3) = y2, say, (14) 


‘ v 2) "in —4) 2 ¢"(n—4) Wn—2) /Gin—2) 2 5 
and OP; Co -—¥# Ce, Ce-® (Cr yz, say. (15) 


Then inversely interpolate from y3—y3 to make this quantity zero (at 
which point y3 = y3 = y?). The relation 


CR? = (y¥s—yCr” (16) 
will give this quantity or its sign if required. 
Sometimes all the last three rows of (2) show loss of significant figures. 
Then accuracy is increased by finding yj, y? satisfying 
Ce-a yi rs Ce-2 yz aa Ce-® = 0, Ce) ys —C@-,9 y+ Ce-» ¥, = 0 
(17) 
and inversly interpolating from y23—y} to make this quantity zero (the 
values of y?, y? close to y2 are those chosen). 
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This loss of significant figures is usually caused by the presence of two 
pairs of conjugate complex zeros with real parts close together and the 
loss is usually more than compensated for by the fact that knowledge 
about two pairs of zeros is being found. In fact, in this case, the alternative 
values of yj, y? in (17) are also found to be close to each other and inverse 
interpolation from their difference will often give a second possible value 
of Xp. 

Bairstow’s (9) method of improving a quadratic factor is, in general, 
a very economical way of proceeding once x, and yz have been found to 
three or four significant figures and it is probably better to use this method 
in the final stages, especially if six or more figures are required. 

When a linear or quadratic factor has been found to sufficient accuracy 
it should be divided out of f(z) and the reduced polynomial factorized. 

Finally, it is perhaps worth pointing out the advantage of the present 
method over Frazer and Duncan’s (2) methods. In their methods the 
test function C{">" (or this multiplied by some quantity) alone was used 
to find a true value of x). Unfortunately C(">” is a rational function of x 
and can change sign through an infinity. Thus, until a good first approxi- 
mation making C\">"” small has been obtained, the choice of the next value 
of x is uncertain and convergence cannot be guaranteed. 


2.3. Approximate times of solution 


If all zeros are complex and have real parts of the same sign then, to 
six-figure accuracy, a competent computer with a modern desk machine 
should be able to solve a quartic in under three-quarters of an hour, and 
a sextic in under two hours. Eighth or higher degree equations are best 
solved by the a method. 


2.4. Example 


We will now solve by the 2 method, the quartic 
z4+-0-93523+- 1-2802z?+-0-425z+-0-265 = 0. 


Now the case x = 0 for which p{0} = 0 shows that one a, must lie in the 
range [0, —0-935/4] = [0, —0-233...], and this is a smaller range than 
[0, —[e,/¢o}!*] 

The value x = —0-1 is conveniently near the middle of this range. 
Again, p{—0-1} = 0 but the value of y2—y? is small. Hence, linear inverse 
interpolation can probably be used and this suggests = —0-126. In fact, 
p{—0-126} = 2 but y3—y} is much smaller and the remaining values of x 
were chosen by linear inverse interpolation. The results are summarized 
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as follows (CY is given in some cases; its value is not important but its 


sign 1s): 
x 
o"9g350 
ol 05350 
0-126 0°4310 
O-1205 045300 


o°12090160 0°451330 


Finally, inverse 
: —0-120923, 


ey 2 
=ly 21+ Yo 


ce 
° 8255 
06088 
06994 


0°697441 
06097550 


% 
0°4545 
03607 
O° 322 
0° 331678 
0331005 


0°3210 
O° 335 
03291 

0° 331136 
0° 330996 


vs—¥a 
O°1335 
o-o251 

— 0°:0067 
0000542 
0°00000G 


c} pix} 
o-r1102 ° 
O-0175 ° 

— 09-0029 2 
T ° 
° 


linear interpolation from the last two results gives 
y2 = 0-330994. The corresponding quadratic factor is 
2? +-0-2418462+-0°345616. 


Divided into the quartic it gives as quotient, z*+0-6931542+-0-776747, 
with zero remainder to six decimal place accuracy. 


The stage x 


)1-0000 (ec, 


(Cy 
(C2 
(Cy )1-0000 (Cys 
(ys =)0-3224 
(ys =)0-3291 


0-0067 


(Cig ° 


)0-9350 
)0-8090 


)0-6830 
)0-5570 


)0-4310 
(= C,) 


(cy =)1-2800 
(cp, =)1-1781 


(cp. =)1-0920 
(¢y3 =)1-0218 
(= C2) 


0-6994 
( = CP) 


(c, =)0-4250 
(C3. = )0-2766 


(¢3.2 =)0-1390 
(= Cs) 


0-1390 
(= Cy) 


— 0-0029 


~ 0-126 is now given in full to show a convenient format: 


(ce, =)0-2650 
(cy, =)0-2301 
(= C,) 


0-2301 
(= CY) 


0-2301 
( - C?) 


(= C}) 

3. The a method 

3.1. Description of the method 

This method consists essentially in factorizing f(z) (see (1)) into quadratic 
factors of the form z*—az+-6, the computer choosing a numerical value 
of a then obtaining two equations for b and continuing until a value a, 
has been found such that these two equations have a common root by; 
z*—a,z+6, is then a factor of f(z). We shall, however, deal only with the 
case in which n is even (= 2m, say), for if f(z) is of odd degree a real linear 
factor can quickly be found and divided out. 

Now, when f(z) is divided by z?—az+-6 there is a remainder of the form 


R,(b)z+ R,(d), (18) 


where the R;, which are polynomials in 6 whose coefficients are polynomials 
in a, are given by 
Rob) = (—)™ [C9 0" — Tony OP 2+ mn —2 O* 2 —.... + 
+(—)™ ano O+(—) Com] (19) 


R,(b) = (—)™ ry yy 0? — 1g. OF +15, mq —p OF —... + (—)™ ama] (20) 








$$ 
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and the r;; are generated from the array 


(Cy N10 (Co =)roo (C3 "3,0 - + + (Com =)Fom—20 (Cam—1 em—1,0 
ria Ton 3,1 Tem—2,1 Tom—1,1 
Ti e2 "3,2 
T1.m-2 T2,.m- 2 T3.m 2 Tam 2 Tsim 2 
Tim-1 Te.m—1 T3.m—1 
Tim (21) 


according to the rule, 


Tu = Tehu-1 TN -1,u (22) 
with c, = 7,,,- This result can be proved by induction without difficulty. 
In words the rule is simply that each r;; is the sum of the term above and 
the product of a and the term to the left. In fact the r;; are identical 
with the c;,; of (9) if a x, but there are fewer (about half) of them for 
the same value of n. It can be seen that if a test using a value of 2 has 
been carried out, then from the resulting c¢,; it is possible to obtain the 
values of r,;, for a having the same value as x, without further work; and in 
fact this combination of the x and a methods is sometimes useful. 

Accuracy can be checked by using one of the relations 


f(l) = R,(a—1)+ R,(a—1), f(—1) — R,(—a—1)+ R,(—a—1) 
(23) 
which can be seen to hold from the fact that 22—az--b is zero if z l, 
b a—l, or2z lL. b a—l. 


When the r;; have been computed the h.c.f. of the two polynomials (19) 
and (20) is found. This can be done by applying Routh’s test to the C; 
where 


’ _ y 2 y nal - Y — y 2 
Cy =t, Tim> C2 =Tem-a C3 = Tam Ce = Tam-2 
y —_ . ’ . » y —. y - » 
( 5 rs.m 2° sf in-3 Tom -2,19 am-1 ~ Tom ~1,1° ( 2m Com: (24) 


The aim is to find a = ay, say, such that the C\*>" of (2) is zero, when, 
simultaneously, 


Wn—2)/¢"(n—2) — "n—3)/¢"(n—3) 30 t7 on 
n ( 2 ( n—1 ( n-3 bo, say. ) 


n -~) 


Further (see section 4) it can then be said that R,(b,) = R,(b,) = 0 and 
hence, that 


~2 


—A,z+b, 


is a factor of f(z). In fact a, may be found in a systematic manner by 
noting the number of changes of sign p{a] say, of the C’” in the Routh test 
array (2) derived from the C; of (24). Once a, and a, have been found such 
that p[a,] + p{a,] then (see section 4) a true value of a, must lie in [a,, a9]. 
Further it can be shown (see section 4) that p[—oo] = 2m and p[ +00] = 0, 
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so that ranges |a,,a,| exist. It does not seem possible to specify such a 
range initially, in general, except in the important cases when p[{0] = 0 
or 2m. In either case it can be proved (section 4) that p[ —c,/(mcg)| # p[O]. 
In other cases one or two trials will usually establish such a range, or by 
one or two applications of the 2 method, a change of origin will be found 
which will convert f(z) to a polynomial all of whose zeros have negative 
real parts. 

As in the x method, there may be loss of significant figures in C{*>" and 
then the last few rows of (2) can be calculated as follows: Evaluate 


b, Ce-* Ce, ce = Ce —b, Ce-, (26) 
b, = C#-*/Ce-2, Ca—) = (b,—b,)Ce-. (27) 


Inverse interpolation from 6,—b6, can be used to find the value of a for 
which b, = b, = by. 

It should be noted that if p[a,] = p[a,| there may nevertheless be true 
values of a, in [a,,a,] and that this can occur if there are more than two 
real zeros of f(z). In this respect the a method differs from the 2 method. 

The number of calculations per stage in the a method is lower, however, 
and it fits in with Bairstow’s method (9) somewhat more easily (this is 
certainly to be recommended once four figures have been obtained if 
n = 8 or more). Finally, though it might appear that real zeros can only 
be found in the form of quadratic factors by the a method, real zeros can 
be found in the usual manner if, in computing the r;;, we perform the one 
additional calculation 


f(@ = Com tOTom—1,1- (28) 


3.2. Approximate times of solution 

If all zeros are complex and have real parts of the same sign then to 
six-figure accuracy a quartic can usually be solved in half an hour, a sextic 
in an hour and a quarter and an eighth degree in two hours and a quarter. 


Example. We now solve the following sextic equation 
2843-3125 +-7-71z4+ 10-9423 + 10-942?+-6-70z+- 2-68 = 0. 


Here p[0] = 0; hence, an ay lies in [0, —3-31/3 = —1-10...]. We choose 
for our second value of a, —0-6 (—0-55 would be more logical but 0-6 
simplified the arithmetic slightly), and 6,—6, is found to be much smaller 
than for a = 0 though pia] is still zero. The next trial value of a is there- 
fore —0-7 and pla] is 2. The value a = —0-64 seems plausible and from 
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this stage onwards inverse Lagrangian interpolation was used. The follow- 
ing is a summary of the main results: 


2 ) 06 —O'7 0°64 — 0°645600 0°64556200 
, 3°3100 1°51 1°2100 1*3g00 1*373200 1°37331400 
Ct?) 4°4050 1°9047 1°5382 1°7597 1*739160 1°73929883 
3 I*1O51 09894 o811g9 0°9438 0°935856 0°93591222 
De 06621 )392 O*9100 0°9 368 0°935904 O'935Q9I1II0 
bs Dg 0°4430 2 O°0gdI 00070 — 0000045 O*OOOOOI112 
pia} ° °) 2 ° 2 ° 
Interpolation finally gives ag = —0-64556286, b, = 0-93591094 and divi- 


sion of the sextic by z?—a,z+-b, gives 

244. 2-664437 1423 +- 5-0560274022 + 5-18363175z-+- 2°86352033 
with remainder —0-00000004z—0-00000000. The stage a = —0-7 is now 
set out in full to show the format. 


(cy = )1-0000 (ec, 33100 (ce, =)7-7100 10-9400 10-9400 6-7000 2-6800 
(T15 )2-6100 (ry, )5-8830 6-8219 6-1647 2-3847 

(Cy =)1-0000  (r,, =)1-9100 (r,., =)4:5460 3-6397 6-1647 2-6800 
Cy /C; 0-8264 (7r,, =)1-:2100 3°6397 2-3847 

{ C*") 
C,/CS = 0-7866 5382 4-1940 2-6800 
(= Of) (= CY) 

b,( 0-2766/0-3407) 0-8119, 0-3407 0-2766 

b.(= 2-6800/2-9451) = 0-9100, 2-9451 2-6800 

b,—b, 0-0981. 


The a method can now be applied to the quartic quotient. We know 
that an a, (of this quartic) must lie in 0, —2-6644.../2 since all real parts 
of zeros are negative. Further, the work on the sextic suggests that a, is 


(algebraically) less than —0-7. Hence a = —1-0 seems a convenient value 
to try and p|—1-0] is found to be 2. Hence a = —0-84 is next tried (6, 
and 5, are very close here), then a 0-83, and from this stage onward 


inverse interpolation can be used. 

After a, and 6, have been found to about four-figure accuracy, say, by 
the application of Lagrangian inverse interpolation using three or four 
values of a near the suspected value of a), Bairstow’s method of 
improving a quadratic factor can be used. Indeed if the computer is 
familiar with these techniques, their combination is a very efficient way 
of using the a method. 


4. Properties of, and proof of Routh’s test (and the fundamental 
theorem of algebra) 
4.1. Examination of the causes of change in p{x} and the accuracy needed 
in its calculation. 


Now the C;, the coefficients of powers of z—z in the expression (8) for 
f(z), can be seen to be polynomials in x and therefore always finite. The 
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C&,, of the Routh test array (2) calculated from the C; are, however, 
rational functions with C=" in their denominators, and hence may have 


infinities. We now consider what happens when C\” goes through a zero. 
Applying the relation (3) we see that 


"r+1) (r—1) se mr) Yr) 9» 

( P+ 1+2s CY +1+238— -{C jCP +28+2: (29) 
r+ (r) ) wn /gr) | “r)\2. q 

( Mu se = CF 24 w—C3 " +oe+aly 4 [Cv} 3+ O(C, ) ’ (30) 
“r+3) “r+1) f bh re) yw « 
Cristes = Cristes—tC A! JOP Sat ae: (31) 


r+1+ 
infinity, but that the C7? ,, are finite and so are the C\"7,°) ,,; for when we 


substitute for the C'’'4)) ., in (31) we find the terms in 1/C%” cancel. Hence, 
since rows of (2) are ound from the two rows above, all subsequent 


Now we see that as CY” goes through a zero, the C'’"5),, go through an 


rows in (2) are finite, in general, unless a higher C\” is simultaneously zero 
with CY. 

Now, (29) shows that as C’” goes from —0 to +0, C7,” goes from + to 
—CO°>PC®, © and C7,” does not change sign. Considering all possible 
combinations of sign of Co», C"), we find that the numbers of changes of 
sign in the sequence C("5"”, Cw”, C7”, Ci? does not alter as CW” goes 
through a zero. Therefore p{x} cannot alter, unless r= »—1 or n (in 
which cases C\"7 is non-realizable). 

If C!) = C, = f(x) goes through a zero then consideration of the possible 
change in signs shows that p{x} can change by one only (if, simultaneously, 
C\"> is zero then a change of two is possible, f(z) then having two co- 
incident real zeros). If C\">" = 0, and simultaneously, C\” is non-zero and 
C(n-2 /O\™—» is positive, consideration of the possible signs shows that pix} 
changes by two as C\"> goes through its zero but does not change if the 
ratio C\"-2) /O\ >?) is negative. Neither will p{z} change if C\">" changes 
sign through an infinity, for then C\"-,” must go through a zero. 

A consequence of the above malian is that a small error in calculating 

\outh’s test will not affect p{xz} even if, say, C\” is calculated to be +e 
when it should be —», both being small, unless r = n—1 or n. In such 
cases x is near the real part of a zero and greater accuracy is necessarily 
required. 

Since the Routh test array (2) is derived from a process for finding the 
h.c.f. of two polynomials, i.e. the odd and even powers of f(z) it follows 
that if, for x = x), C\">” is zero, then all polynomials 


Cr) y— —ce;” 2y"~ 24. C= r) ay’ *—... (32) 


have some common factor, y?—y3, say. In particular 


Cin—2 y2—Cn-® = 0, Ce-) y8—Ce-—» y, = 0. (33) 


aa 
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Also, putting r = 0 and 1 in (32), we see that 
Cy ya—C,yn—? | CC. —... == GU, C,yg~'—C,y3-* +... =2 ©. (34) 


ence, if z—1% = +ty, 8) and (34) show that f(z) is zero. Thus 
H f ° -iyYo, then (8) and (34) show that . ro. Thu 
a change in sign through a zero of C\">), and, hence a change of two in 


rT 
f 


pix}, is shown to be coincident with the existence of a pair of complex 
conjugate zeros. 

Similar observations apply to the C,; of (24), with y? replaced by 6. 
Here, as a goes through a), C\">” goes through a zero, pia| changes by 
two and R,(b) and R,(b) have a common factor b—b,. Hence from (18) 
it follows that f(z) divided by z?—a,z+, gives a zero remainder. 


4.2. Negative terms in the Routh’s test array 

If all zeros have negative real parts f(z) is generated by factors of the 
form (z+-A) and z*+- Bz+-C where A, B, C are necessarily positive. Hence, 
if cy is positive in f(z) (see (1)) so are all the c;. Conversely if one c; is 
negative then p|0} ~ 0, and at least one of the C’” must be negative also. 
Now each row of the array (2) is generated from the two rows immediately 
above, always in the same way. Hence, if in carrying out Routh’s test 
to calculate p{x} or p|a| we find one of the terms in a row to be negative 
then necessarily one of the subsequent C” must be negative. Thus p{zx} 
cannot be zero, and, if this is what was wanted, no further work need be 
done. 


4.3. Proof of Routh’s test (and the fundamental theorem of algebra) 


Now, ifc, = 1, as x > +, the C;, of (8) can be seen to have the asymp- 
totic values 
Uy = ® C; ~ (") 


Applying Routh’s test to the C, we find the terms in the array (2) have 
the general form 


, n 
CP) am , 
aa 2r+ 2s 


[s+ 1)(8+-2)...(84 r)}{(m+-1)(n+3)...(m+ 2r— 1) ]27-1a*r +28 


© [(2r+- 28+ 1)(2r+ 28+-3)...(4r+ 2s—1) ][n(n—2)(n—1)....(n—2r 2))’ 
(35) 





n 
"(2r+1) - 7. 
( 2r+1+2s ( ) x 


2r+2s-+ 


[(s+1)(s+2)...(s4 r)][(n + 2)(n-+ 4)...(m+- 2r) |27-1y2r +2841 





*[(2r+- 28+ -3)(2r+ 28+-5)...(4r+ 284 1) || (n>—1)(n—3)...(m —2r+1)] 
(36) 
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(these satisfy the relation (3) and have the right values when r = 1). 
Since, in any C}", j is less than or equal to n, the coefficients of powers 
of x in (36) and (35) will be seen to be essentially positive. Hence, C” has 
the same sign as x’, and so p{-+-00} = 0 and p{—oo} = n. Hence, by 4.1, it 
follows that there exist at least n zeros of f(z) which are, in general, distinct. 
It is well known that f(z) cannot have more than n zeros and so it has 
exactly n. Then, if p{x,} = p, it follows there are p real parts between z, 
and oo, which proves Routh’s test. 


4.4. The application of Routh’s test to the a method and a further proof 
of the fundamental theorem of algebra 
Here, n = 2m, and the C; of (26) can be seen from (23) to be poly- 
nomials in a, with asymptotic values, as a + 00, 


ao anf aa 6-CSIP 


+6§ 
Cy~ re ‘at Cs ~ P ' Ja’, tte (37) 


The same argument as in 4.1 shows that since the C;, are finite, p[a] can 
change only if C2">" or C2” (calculated by (2) from the C;) go through 
a zero. But C2” = C,,, = Cy, = constant. So change can occur only if 
C2™—>) goes through a zero (and therefore will be by two). 

Now the C'” of (2) have the asymptotic values 


(2r) 


2r+ 2s 


[ (2m —1)(2m—3)...(2m—2r—1)][s(s+1)...(s+-r—1) ]a”+* 
| (2r+-28-+ 1)(2r+- 28+ 3)...(4r+28—1) |[m(m-+ 1)...(m+r—1)] 


m-—--T--s x 
2r+2s+-1 


» (38) 





[(2m+-1)(2m-+3)...(2m+ 2r+ 1)][s(s+1)...(s-+r—1) ja® +2841 


39 
. [(2r-+-28+3)(2r+28+5)...(4r+ 28+ 1) || (m—1)(m—2)...(m—r)] —_ 





and hence the C” will have the same signs as a’ and so p{|+oo] = 0, 
p| —%0] = 2m. Thus at least m factors of f(z) with different a’s exist 
(which incidentally again proves the fundamental theorem of algebra). 

Now, if p[0] = 0, all real parts are negative. But (23) shows that 
C, = ¢,+macy. Hence if a is (algebraically) less than —c,/m¢o, C, will be 
negative, and, hence, p[—c,/mcy] #0. A similar argument holds when 
p[0] = 2m, which shows that p[—c,/mc,] A 2m. 

5092 .47 ce 
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DUAL FOURIER—BESSEL SERIES 


By J. C. COOKE (Farnborough) and 
C. J. TRANTER (Royal Military College of Science, Shrivenham) 


[Received 9 June 1958] 


SUMMARY 
A method is given for determining the coefficients a, in the ‘dual’ Fourier—Bessel 


series - 

- 
yona 
=1 


Ji(a,r) = F(r) (O<r< 1), 


n™y 
n 


x 

¥ a, Jar) 0 (l<r< a), 

n=1 
where — 1 P 1, F(r) is specified and «a, is a positive root of J,(a,a@) = 0. Such 
series arise in the solution of potential problems in which the application of a finite 
Hankel transform is appropriate and in which one of the boundary conditions is 
a ‘mixed’ one. As an illustration, the method is applied to find the potential due 
to an electrified circular disk situated inside an earthed coaxial infinite hollow 
cylinder. This shows that useful results can be obtained provided that the radius 
of the cylinder is sufficiently large. 


1. Introduction 

POTENTIAL problems in a semi-infinite medium in which different condi- 
tions hold over different parts of the same boundary can often be con- 
veniently reduced by the use of a Hankel transform to the solution of 
a pair of dual integral equations. When the medium is confined within 
a circular cylinder or between a pair of parallel planes, such problems 
can be similarly reduced by the use of a finite Hankel transform to the 
determination of the coefficients in what may be termed ‘dual’ Fourier- 
Bessel series. In the physical problems so far considered, these dual series 
are of the form 


x 


> a, 5(a,7) = Fir) (0<r <1) 

n=1 (1) 
> a, J,(a, 7) = 0 (l<r<a) 
n=1 


where —1 < p < 1, F(r) is specified, «, is a positive root of J,(a,a) = 0 
and the coefficients a, are to be determined. 

The purpose of this note is to give a method for determining the coeffi- 
cients a, from equations (1). The solution given is formal only and the 
difficult question of convergence is not considered. The method used is 
similar to that previously used by one of us (1) in the case of dual integra 
equations and, under certain circumstances, the solution given is suitable 

{[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 3, 1959] 
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for computation purposes. As an illustration, we consider the problem 
of an electrified circular dise situated inside an earthed coaxial infinite 
hollow cylinder. 


2. The reduction to a set of algebraical equations 
[t is sufficient for the purpose in view to take v to be real and it follows 
from Tranter (2) that, if m is zero or a positive integer and vy > —1, 


x 


> Jysamsr+ipl%nMl%nt) _g (1 <r <a). 
* ip J? | Xn a) F 


By taking S oJ (a), 


1+tp J2 / m”“y+2m+1+ip 
yt te J? (a, a) a, 


the second of equations (1) is therefore automatically satisfied. 

The coefficients b,, have now to be chosen so that a,, satisfies the first 
of equations (1). Substituting from (2) and interchanging the order of the 
summations, 


S bn, S 


m=0 n=1 


J, +2m+1+ tp(%n J, Xn r) 
yt 4P J? (x, a) 


= F(r) (0<r <1). (3) 


Now (3), if s is a positive integer or zero, 
J) 428+1+4p(%n) 

1+4p 

xt a7 


x 1 
C(v+s+1) % 


* WLI) +1+ yp) 


rv+1(1 —r?)'P A (1+-4p+-y,v+ 1, r?)J,(a, 7) dr, 
0 
where 

Fi(i+4pt+y,vt+1,r?) = .F(—s8,1+4p+v+s; v4 
is Jacobi’s polynomial. Hence multiplying (3) by 

r’+1(1—r?)tPF(1+4p+y, v+1, r?), 

integrating with respect to r between 0, 1, interchanging the order of 
integration and summation and use of (4) gives 


x 


“~ J 40 +144p(% VW) 994144p(% ) 7 5 
U m ? ni” v+2s Dvn = E » 8, ‘ f 
22 J? 4 (%y, 4) stl ei 

where, for brevity, we have written 
E(v, 8, p) 
1 
T(v+s+1) ‘ ‘ 
NS A. a Sl Tee, Se vtl(] _»2)§P F(]4+-1ytyp. y+, r?) Fir) dr. 
T+ 1) (s+1+4p) | a a Ei os al a a fa 
0 


(6) 
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Equation (5) with s = 0, 1, 2, 3,... provides a set of algebraical equations 
to determine the coefficients 6,,. Once the values of b,, h->ve been found 
from this set, the coefficients a, can be evaluated from equation (2). 


3. The iterative solution of the algebraical equations 

The set of algebraical equations (5) can, under certain circumstances, 
be solved by an iterative process as follows. Provided v is not a negative 
integer and that vy > —1—4p,, it is known (Tranter (2)) that, for all values 
of m and s under discussion here,f 

2 ~ J, 2m+1 +ip(%n J, +2s+1+ ip(%n) oo _- i 


2 / 2 2 » i ae 
a® x2 J? (x, @) 2v+-48s+2+p 





@ 


aes K,(t) t t = 
-(—1)™ *=sin$pa | i) Iysamsrsto (Hosa ao(2) dt, (7) 
0 
where 5,,,, = 0 when m #8 and 6,,, = 1 when m = s. Writing 
LS) 
a ae  K,(t) t t 
( ae l | a +8 + 1(2y +4s8+- “+p)-om 4pr | Tip bvsamst stn) loszensto( dt, 
0 
(8) 

and substituting from (7) in (5), the algebraical equations to determine the 
coefficients b,, become 


x 


9 
b,— > by, LEP = —(2v4+48+2+p)H(v,8,p) (s = 0,1,2,3,...). (9) 
m=0 a 
The iterative solution of this set of linear equations (which can be shown 
to converge if a is large enough although it is difficult to determine precise 
limits for this) is given by 
b,= > 6” (s = 0,1,2,3....), (10) 
r=0 
2 
where 6° — a2 (ev + 48+ 2+ p) Ely, 8, p) 
(11) 


and 6”) = > Emin bm» (r = 1,2, 3....) 
m= 


+ The particular case of vy = —}, p = —1 is of some importance in practice. In this 
case the inequality vy > — 1—4p does not hold but the result (7) remains valid except when 
m = 8 = 0. By modifying the analysis given in (2), it can be shown that the result for this 
case is 


2 a2 J? = 
a pores | Xn 3(%n 2) 


x a 
2 T2a) g.(%)—2 Ti(t/a)—1 


t(e** +11) 
0 
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4. The evaluation of LL”). 


m.,8 

The theoretical determination of the coefficients a,, of the dual series (1) 
is contained in equations (2), (10), (11), (8) and in a practical problem, 
the chief difficulty lies in the computation of L&”). By expanding the 
Bessel functions of argument (t/a) in equation (8) in ascending powers of 
t/a, Li? can be expressed as a series of descending powers of a with the 


aid of the integrals 


PK (t) 
AM) ay 
L(t) 


M® = (12) 


These series are quite rapidly convergent if a is fairly large compared with 
unity. 


TABLE 1 


Values of M”? = 





° 





1*36768 


0°64689 


2°06986 


16°168 


231°69 














5286°6 





The values of v occurring most frequently in practice are +4, 0, and 1 
and the integrals M‘°) have been evaluated by Knight (4) for n = 0, 2, 4, 
6, 3, 10. These values, together with some of those for M‘) (which have 
been calculated in a similar manner) are given in Table 1. These integrals 
are laborious to compute and, to make the most of the method, Table 1 
could be usefully extended to higher values of n. Useful approximating 
polynomials for K,(t), [,(t) when v = 0 and 1 have recently been given 
by Allen (5) and these should reduce the labour of extending this table 
should such a project be undertaken. Values for M‘-» and M‘ are also 
given in the table—these integrals can be expressed in terms of the 
Riemann zeta function and are easy to compute from a table of that 
function. 
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Some inequalities involving the integrals M® are useful in estimating 
the error when LY?) is computed by the above method. Knight (4) gives 


me-*(2n+-1)! 
Q2n+1 


M® < 2+log,2—y+ (13) 


where y = 0-577215... is Euler’s constant. For n > 3, a better inequality 


™ 3-15(2n)! 


MY < (14) 


Q2n+1 


It can also be shown that, for n > 1, 


8(2n)! 
(1) 

M an < 92n+1 (15) 
and it is almost certain that this can be improved upon. 


As an example, using equation (8) and Table 1, it can be shown that 


= 0:33818a-3+0-10190a-5+-0-05426a-7 + 0-03457a-%, 
with an error less than =" 
5. An example 
As an illustration, we consider the potential due to an electrified circular 
disk charged to unit potential and situated inside an earthed coaxial 
infinite hollow cylinder. An approximate solution to this problem has 
been given by Smythe (6) working by quite a different method. Using 
cylindrical coordinates (r,z), we take the disk as z = 0, r < 1 and the 
cylinder as r = a, —c <z< oo. The potential V has to satisfy 


@2V 12eV. av 
= — = Q, 16 
or? ' r Or 2? (16) 


with V=1,0<r <i: ov =01<r<a, when z = 0, (17) 
since there is symmetry about the plane z = 0, and 
V=0 when r=a, V>O as jz! >. (18) 
If V is the finite Hankel transform (of order zero) of V, i.e. if 
a 
v= | rVJ(ar) dr, (19) 
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where « is a positive root of the equation J,(aa) = 0, then, by the usual 
procedure (see, for example, Tranter (7)), equations (16) and (18) lead to 
(20) 

The solution of (20), finite at infinite axial distances is 
V = A(aje-**|, (21) 
where A(a) is a constant depending only on «. The inversion theorem 


appropriate to (19) then gives 


(22) 


Writing «, A, = a, J3?(«,a) and substituting from (22) in the ‘mixed’ 
boundary condition (17), the coefficients a,, have to be found from 


n=1 


> Xn Ma, Jy(xp, r) ha? (O 7 : : | 


x 


> a, A(a,7) = 0 <r<a)| 
1 


n 


These are a special case (v = 0, p = —1, F(r) = 4a*) of the more general 
dual series already discussed. 

The physical quantity of most interest in this problem is the electro- 
static capacity of the system and here this is equal to the total charge 
on the disk. The surface density of charge on the disk is (—}7)(éV /éz),_9 
and, taking into account both sides of the disk, the total charge is 


Using equation (22) and the relation between A, and a,, this can be 
written 
. @ 


rJj(a,r)dr, i.e. = 3 Gn Jy(%p)_ 


2 
= n=1 “n 


if we assume it is legitimate to interchange the order of integration and 
summation. Putting vy = 0, p = —1 in (2) 


(a, x, 2) > Fam-+i( %n) 
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and hence the capacity C, assuming an interchange in the order of the 
summations is permissible, is given by 


ct : Toy 4(%n (Xp) 24 
>, bm 4 at J3(a,a) (24) 


Now it is known that, Tranter (2), 
a 


-> Fam4(%nJo(%n?) __ [ thJon4(t)Jo(rt) dt, 


2 Z 2 
a hot, xt JF( xX, a) J 


0 
so that, multiplication by r and integration with respect to r between 0 
and | gives 


Pe Jom, + 3( Xn J (ay) 


an } t-¥Jom, «(tdi (t) dt 
é 


From Watson (8) the infinite integral on the right vanishes except when 
m = 0 and then its value is (2/7). Equation (24) then shows that the 
capacity is given by 2\4 
Cn (=) by: (25) 
\7T 


With v = 0, p = —1, F(r) = 4a*, equation (6) gives 


a* T(s+1) 


E(0,8, —1) = 3 Fed) 


1 
| r(1—r?)-*.F,(4, 1, r?) dr. 
0 
Now A,(4, 1,r?) = 1 and the orthogonality relation for the Jacobi poly- 
nomials (9) shows that 

a2 
(27)¥” 


and the first of equations (11) then gives 


E(0,0,—1) = E(0,s,—1) = 90 (8s > 0), 


7 


9 
b(o) ss (=)! 59 =O (es = 1,2, 3....). 


By expanding the Bessel functions of argument (t/a) in equation (8), it is 
apparent that L{\°>” is of order a~®”-**-! and, if we reject terms above a~, 
equations (11) and (25) yield 


o = (=)fa, = 20-159). (26) 
7 7 


Working as in section 4, it can be shown that 


a 2\"[1-36768 | 0-21563  0-09199 
Lis? = (-) |r ee 





7 


and these relations lead to values of the capacity of 2-814/7 when a = 2 
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and 2-3255/7 when a = 4. These results may be compared with Smythe’ 


estimates that the capacity lies between 2-807/7 and 2-818/7 when a = 2 
and between 2-3253/m7 and 2-3256/7 when a = 4. 


“I 
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DRAG FORCE ON A GOLF BALL IN FLIGHT AND ITS 
PRACTICAL SIGNIFICANCE 


By D. WILLIAMS (Royal Aircraft Establishment, Farnborough) 
(Received 25 April 1958] 


SUMMARY 

An experimental method of measuring the aerodynamic drag of a golf ball in free 
flight under typical playing conditions is described and the results obtained are 
discussed. The method depends on deriving a curve in which the ‘carry’ of the 
ball is plotted against initial ball velocity off the tee, so that only these two para- 
meters need to be measured. A main conclusion is that the variation of the drag 
coefficient with Reynolds number is such as to make the drag vary linearly rather 
than as the square of the velocity (over the relevant speeds), with a consequent 
advantage to long hitters. 


1. Introduction 


Tue forces acting on a golf ball in flight have been the subject of a number 
of investigations in the past. Of these probably the latest is that carried 
out by J. M. Davies (1) in 1949. Apart from presenting a comprehensive 
review of previous work, he describes in his paper a number of experi- 


ments carried out in the small wind tunnel of the B. F. Goodrich Co. 
The effect of various degrees of surface smoothness and of rotational 
velocity on the lift and drag of the ball were measured and some interesting 
results were obtained. Unfortunately, however, all the experiments were 
made at a constant tunnel speed of 105 ft/sec and therefore at a constant 
Reynolds number of 88,000. It was therefore not possible to estimate 
how the aerodynamic forces on the ball vary with the higher Reynolds 
numbers achieved in practice. 

It is true, of course, that a fair amount of data on the drag of spheres 
of varying degrees of roughness over a wide range of Reynolds number 
already exists and is recorded in the standard textbooks. It is all based, 
however, on wind-tunnel tests on non-rotating spheres with a type of 
surface roughness very different from that of the conventional dimpled 
golf ball. But, although it would be hazardous to use this data for estimat- 
ing the forces on a golf ball in free flight, it has amply demonstrated the 
interesting fact—and particularly interesting from the golfer’s point of 
view—of the sudden and very marked drop in the drag coefficient at a 
certain critical Reynolds number. 

It is now well known that this drop is due to the fact that the transition 
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from laminar to turbulent flow occurs earlier in the passage of the air over 
the ball, thereby narrowing the wake and reducing the drag. 

Past experimental work has also shown that the smoother the surface 
of the sphere the more sudden the drop in the drag coefficient, and that, 
beyond a certain degree of roughness, the drop starts at a much lower 
speed and is distributed over a much wider speed range. Again it is 
difficult to translate the results of experiments on non-rotating spheres 
with particular types of roughness, in tunnels of a particular degree of 
turbulence, in terms of the spinning golf bail in free flight with its own 
peculiar type of surface roughness. 

The ordinary handicap golfer’s interest in this problem arises from his 
suspicion that the long drives achieved by the professional player are 
largely due to some peculiarity in the behaviour of the ball which he 
himself is unable to exploit but which he feels is ‘just round the corner’ 
if only he could achieve that little extra initial velocity. The matter is 
also of interest to the ruling body that governs the game, for if it were 
found that there is a critical speed which, if reached, gives the player some- 
thing for nothing, so to speak, some action might be taken to obviate 
what might be considered an unfair advantage. 

It was with such considerations in mind that the writer evolved a fairly 
reliable method of determining how the drag of a golf ball varies with its 
speed under typical playing conditions. To make the conditions typical, 
the ball must be struck by a lofted club-head similar to that of a driver 
so that, with typical back-spin, it follows the orthodox flight path of a 
gor xd drive. 

Two parameters alone need be measured—the initial velocity of the ball 
and its ‘carry’, i.e. the distance it covers before striking the ground. By 
plotting the ‘carry’ against the initial velocity a curve is obtained the 
tangent to which at any particular initial velocity gives the aerodynamic 
drag at that velocity. 

It is claimed that this is not a ‘rough-and-ready’ method of obtaining 
the drag coefficient of a golf ball under typical playing conditions but 
a perfectly sound method that leads to results that are inaccurate only 
in so far as the two relevant parameters are inaccurately measured. This 
is the only reason for the use of the phrase ‘fairly reliable’. 

The necessary tests, which covered a range of carry extending from 150 
to 230 yd, were carried out by Messrs. Dunlops (through the courtesy 
of Mr. S. G. Ball), and the data they obtained is plotted in Fig. 1. From 
this it is seen that, between the chosen limits of initial velocity, the curve 
is well represented by a straight line, so that the tangent to the curve is 
constant. 
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2. The theoretical basis 

Fig. 2 shows the trajectory for a typical drive, in which the initial angle 
a of the flight path is very shallow and its final angle 8, where it meets 
the ground is fairly steep. 

Let 2, be the total carry, v the ball velocity, and v, the initial velocity 
of ball. 

Consider the trajectory AB for which the initial velocity is v, and the 
total carry 2%). Suppose now the ball is teed up at A’ a very short distance 
5,2) behind A and slightly below it by an amount 5, x, tana. If, starting 
from A’, the ball is given an added velocity 5v, that sends it along precisely 
the same path AH B as before, it will reach a point B’ on the same level 
as A’ and beyond B by the distance BC’, where 


BC’ = 5,% = B'C’ cotB = 8,2, tan acot P. (1) 
The total carry on the level is thus 
A’B’ = (a +82) = (ao +8, %9 +522) (2) 


and the increase in carry is 


b2y = 5, %) +5, 2% = 5, 29(1-+-tan a cot). (3) 
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Since, in order to follow the same trajectory, the velocity v, at A must 
be the same in both cases, the velocity (vy+4v,)) at A’ must have dropped 
by the amount 5v, in its passage from A’ to A. We can therefore write 


bv dv 
= —— (4) 
01% ox 
or, by (3), and introducing differentials, 
dv dv 
(1+tan acot 8)—° = ——. (5) 
dx, dx 





Now x is normally about 6° and f about 45°, which makes 
(1+tan a cot 8) 1-105, (6) 


and it is clear that small errors in « or 8 are of minor importance. 
Reverting now to the experimental results of Fig. 1, we see that the 
constant slope as measured is 
dx,/dv, 120 x 3/80 4-5 sec, (7) 
and therefore, by (5) and (6), 
(dv/dz), (dv,/dx,) 


. Pe —s -I9AR geac-l 
measured ~ 1-105 = — 5 = —(-246 sec. 


te 


(3) 

The equation of motion of the ball in its line of travel as it leaves the 

tee (if we neglect gravity and the lift due to the spin, both negligible at 
this stage and also tending to cancel each other) takes the form 


dv . 
mv— cos a+ Kv? 0 (9) 
dx 


or, with sufficient accuracy, since cosa ~ 1 


do, Ke _ 
dz'm 


0. (10) 
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Here v = velocity of ball in ft/sec, 

m == mass of ball of weight 1-62 oz = (1-62/16)/g slugs (g = 32-2), 
K = 4CppS = 17x 10-Cp, 
Cp = aerodynamic drag coefficient, 

S = projected area of ball = }nd* = 0-0143 ft?, 

d = diameter of ball = 1-62/12 ft, 

p = mass of air per cu ft = 0-00237 slugs (at sea-level). 


Inserting the appropriate numerical constants in (10) we obtain the 
relation 


Uae 
v\dx 


D 


= —0-00536C,, (11) 


0-246 46-0 
~ (0-00536)v oe 


or, from (8) (12) 

Thus, for ball velocities ranging from 150 ft/sec to 230 ft/see, which 
represents the practical range, the drag coefficient varies inversely with 
the velocity. The total drag therefore increases only linearly with the 
velocity instead of as the velocity squared, as it would for a constant drag 
coefficient. The penalty suffered by the hard hitter is therefore substan- 
tially reduced. 

Using the value of Cp given by (12) in the expression for K, we obtain 
the total drag on the ball as 


(17 x 10-*)C,, v? Ib = 0-000783v Ib. (13) 


Thus at the initial velocity of 190 ft/sec, which gives a carry of 180 yd— 
a typical distance for the average player—the drag coefficient from formula 
(12) has the value 0-242. For a long drive of 232 yd carry the initial 
velocity is 225 ft/sec at which the drag coefficient is 0-204, a reduction of 
16 per cent. 

It is of interest to see how much of the difference in length between the 
long hitter and the average player is due to this drop in the drag coefficient. 
To determine this we use the numerical values just quoted and suppose 
that the drag coefficient for speeds between 190 ft/sec and 225 ft/sec 
remains constant at 0-242—the value appropriate to 190 ft/sec. From (11), 
on putting 0-242 for Cp, we have 


1 dv _ _ 0.005360, = —0-0013 (11) 
vdx 


for which the solution is 


log,v = —0-00132-+-constant. 
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Using the condition that v = v, when x = 0, we write 
log, U/U = —0-0013z. 


The distance travelled while the velocity drops from v, (= 225 ft/sec) to 
v, (= 190 ft/sec) is therefore 


225 
; - 4 = 129 ft = 43 yd. 
(<<oo13)!& rao 129 ft 43 yd (14) 


Comparing this with the actual distance of (232—180) = 52 yd, we see 
that the long hitter gains just over 9 yd. It can be said therefore that 
only some 17 per cent of the long hitter’s advantage in distance can be 
attributed to the drop in drag coefficient. The rest apparently he gains 
by honest effort. 
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